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ABSTRACT 


An  experimental  investigation  was  conducted  to  develop  a  method  of  predicting 
cylinder  indicated  torques  in  a  reciprocating  engine  by  measurement  of  crankshaft  angular 
velocity  fluctuations.  Cylinder  indicated  pressures  were  measured  for  all  three  cylinders 
of  a  two-stroke  Diesel  engine  with  pressure  transducers.  Time-resolved  angular  position 
was  measured  at  the  crankshaft  front  and  at  the  flywheel.  A  six  degree-of-freedom 
torsional  crankshaft  model  was  developed.  Two  solution  methods  are  described  to  solve 
the  equations  of  motion:  a  time-marching  ODE  solver,  and  a  Finite  Element  solution  in 
the  time  domain.  Using  these  methods  with  the  measured  cylinder  torques,  the  angular 
positions  are  predicted  and  compared  to  measured  angular  positions  for  model 
calibration.  An  inverse  solution  method  was  developed  to  determine  the  cylinder 
indicated  torques  from  the  measured  angular  position  at  the  crankshaft  endpoints.  The 
method  is  theoretically  demonstrated  to  be  useful  for  explicit  solutions  for  two-stroke 
engines  up  to  three  cylinders,  and  four-stroke  engines  up  to  four  cylinders.  Experimental 
results  show  that  the  method  is  useful  in  predicting  cylinder  indicated  torques  from 
angular  velocity  measurements. 
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I.  INTRODUCTION 


A.  MOTIVATION 

Reciprocating  engines  require  a  maintenance  system  to  ensure  readiness 
throughout  operational  life.  Currently,  particularly  in  military  uses,  this  is  accomplished 
with  a  regularly  scheduled  maintenance  (RSM)  system,  where  parts  are  checked,  and 
replaced  if  necessary,  and  components  are  serviced  at  intervals  based  on  historical 
expected  failure  rates.  Therefore,  maintenance  is  normally  performed  well  before  it  is 
actually  necessary,  because  the  true  condition  of  the  engine  is  unknown.  Significant 
savings  could  be  realized  with  the  use  of  a  condition  based  maintenance  (CBM)  system. 
Such  a  system  requires  a  means  to  monitor  engine  health  during  operation,  or  with 
operational  tests  that  do  not  require  significant  work  on  the  engine. 

Several  classes  of  faults  that  occur  in  Diesel  engines  can  be  detected  and  localized 
by  measurement  of  individual  cylinder  firing  pressures.  Examples  include  loss  of 
compression  ratio  due  to  cylinder  leaks,  and  improper  combustion  of  fuel  due  to  injector 
problems.  Monitoring  cylinder  firing  pressure  is  an  excellent  means  of  condition  based 
maintenance.  However,  due  to  the  harsh  environment  in  the  cylinders,  the  pressure 
transducers  required  are  very  expensive  and  short-lived.  While  direct  measurement  of 
cylinder  pressures  for  performance  monitoring  is  feasible  and  sometimes  used  in  an 
operational  engine,  it  is  expensive.  Of  course,  a  possible  solution  to  this  problem  would 
be  the  development  of  cheap,  reliable  pressure  transducers  for  use  in  operational  engines. 
Barring  this,  an  alternative  solution  is  the  use  of  indirect  methods  for  estimating  cylinder 
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pressures,  such  as  the  measurement  of  crankshaft  angular  velocity  fluctuations  along  with 
an  appropriate  scheme  for  inferring  the  pressure  waveform. 

B.  STATE  OF  THE  ART 

The  angular  velocity  of  a  reciprocating  engine  contains  small  fluctuations  due  to 
the  variations  of  cylinder  pressures.  In  general,  the  engine  speeds  up  after  a  cylinder 
fires,  then  slows  down  as  the  next  cylinder  is  compressing  in  preparation  for  combustion. 
The  flywheel  is  intended  to  reduce  the  magnitude  of  these  oscillations,  but  they  are  still 
present  and  represent  a  speed  variation  of  several  percent,  which  is  a  measurable  amount. 
A  number  of  researchers  have  investigated  the  possibility  of  predicting  the  cylinder 
pressure  variation  by  measurement  of  these  small  speed  oscillations. 

Freestone  and  Jenkins  [Ref  1]  measured  crankshaft  velocity  with  a  proximeter 
mounted  at  the  flywheel  ring  gear  teeth.  They  developed  a  lumped  crankshaft  model, 
using  inertial  torque  to  account  for  the  reciprocating  piston  masses.  This  model  was  used 
to  develop  a  calculation  of  the  total  gas  torque  in  the  engine  cylinders  as  a  function  of 
crank  angle.  Noting  abnormally  low  peaks  in  the  pressure  waveform  and  their 
corresponding  crank  angle  localized  faults  in  individual  cylinders. 

Mauer  and  Watts  [Ref  2]  measured  angular  velocity  at  both  ends  of  the  crankshaft 
by  placing  a  proximeter  at  the  flywheel  ring  gear  teeth  and  at  a  corresponding  gear 
mounted  on  the  pulley.  The  phase  difference  between  the  two  encoders  corresponded  to 
an  instantaneous  measurement  of  the  total  crankshaft  twist,  which  was  considered 
proportional  to  crankshaft  torque.  No  mathematical  model  was  used,  so  detection  of 
faults  was  realized  by  comparing  the  measured  signal  to  a  signal  recorded  on  a  healthy 
engine.  As  expected,  the  twist  signal  had  peaks  that  varied  depending  on  which  cylinder 
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was  currently  firing,  because  cylinders  farther  from  the  flywheel  cause  a  greater  twist  for 
the  same  combustion  pressure.  Mauer  continued  this  work  [Ref  3],  developing  a  lumped- 
parameter  engine  model  to  isolate  cylinder  specific  torque.  In  this  method,  he  computed 
the  cylinder  specific  torque,  defined  as  the  integration  of  the  total  crankshaft  torque  from 
TDC  of  the  cylinder  in  question  to  TDC  of  the  next  firing  cylinder.  In  a  healthy  engine, 
cylinder  torques  will  all  be  close  to  the  mean,  but  a  faulty  cylinder  will  show  a  reduced 
specific  torque  compared  to  the  other  cylinders. 

Citron,  et  al.  [Ref  4]  used  a  four  degree-of-freedom  model  of  the  engine-drivetrain 
system  that  differentiated  between  the  flywheel  and  the  engine.  Reciprocating  masses 
were  accounted  for  by  an  inertial  torque  component  and  the  crankshaft  speed  was 
measured  with  a  proximeter  at  the  flywheel  ring  gear  teeth.  Total  cylinder  gas  pressures 
were  reconstructed  by  solving  the  equations  of  motion.  Individual  cylinder  pressures 
were  inferred  by  assuming  that  the  majority  of  the  net  torque  at  a  particular  point  was  due 
to  the  cylinder  undergoing  the  power  stroke. 

Connolly  and  Yagle  [Ref  5]  used  a  lumped  engine  model,  assuming  the  total 
inertia  of  the  crankshaft  components  as  a  single  mass.  An  inertial  torque  accounted  for 
reciprocating  masses  and  the  angular  velocity  was  measured  with  a  proximeter  at  the 
flywheel  ring  gear.  A  nonlinear  differential  equation  relating  combustion  pressure  to 
angular  velocity,  derived  from  the  torque  balance  equation,  was  reformulated  to  a  linear 
first-order  differential  equation  relating  pressure  to  the  square  of  the  angular  velocity. 
Connolly  revisited  the  issue  [Ref  6]  to  reconstruct  cyclic  pressure  variability  from  the 
crankshaft  angular  velocity. 
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Lim  et  al.  [Ref  7]  predicted  cylinder  pressures  in  a  four-cylinder  four-stroke  spark 
ignition  engine  by  making  several  assumptions  about  the  cylinder  pressure  as  a  function 
of  the  crank  angle.  The  intake  and  exhaust  strokes  were  assumed  reference  values  (intake 
manifold  pressure  and  exhaust  backpressure)  and  the  compression  stroke  was  estimated 
as  a  polytropic  process  for  each  cylinder.  From  knowledge  of  the  crank  angle  and  the 
firing  order,  power  stroke  pressures  were  estimated  for  each  cylinder  from  the  measured 
angular  velocity  and  the  known  load  torque.  The  method  implied  a  lumped  crankshaft 
model. 

Iida  et  al.  [Ref  8]  measured  angular  velocity  with  a  proximeter  at  the  flywheel, 
and  included  a  correction  for  tooth-to-tooth  variation,  determined  by  measurement  of  the 
flywheel  at  a  constant  rotational  speed.  A  lumped  engine  model  was  used,  in  which  an 
equation  related  the  total  engine  inertia  and  rotational  acceleration  to  the  composite 
torque  applied  to  the  crankshaft.  Integration  of  this  equation  over  a  cycle  yielded  a 
relation  to  determine  Indicated  Mean  Effective  Pressure  (IMEP)  from  the  total  engine 
inertia  and  the  square  of  the  change  in  the  angular  velocity. 

Taraza  [Ref  9]  developed  a  linear  high-fidelity  model  of  a  multicylinder  engine. 
Using  this  model,  angular  velocity  and  angular  deflection  were  predicted  for  the  front  of 
the  crankshaft  and  compared  to  measured  values  for  an  inline,  four-stroke,  four-cylinder 
Diesel  engine,  in  order  to  verify  the  model  parameters.  He  then  determined  harmonic 
orders  for  the  crankshaft  model,  and  conducted  experimental  measurements  at  these 
speeds.  His  results  show  good  agreement  between  the  measured  and  predicted  harmonic 
order  amplitudes.  His  conclusion  was  that  the  measured  amplitudes  of  certain  harmonic 
orders  of  angular  motion  could  be  used  to  determine  engine  mean  indicated  pressure. 
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Additionally,  his  work  demonstrated  the  usefulness  of  a  high-fidelity  crankshaft  torsional 
model. 

Additional  work  on  this  subject  [Refs  10-16]  generally  used  a  lumped  crankshaft 
model  with  angular  velocity  measured  at  one  point. 

Previous  work  at  NPS  by  Bell  [Ref  17]  and  Hudson  [Ref  18],  on  the  same  engine 
used  in  this  study,  demonstrated  that  there  is  information  present  in  the  crankshaft 
rotational  speed  of  a  reciprocating  engine.  Hudson  developed  a  high-fidelity  model  of 
the  engine  crankshaft,  then  used  measured  pressures  to  predict  speed  fluctuations  for 
comparison  with  actual  speed  fluctuations  at  the  crankshaft  nose.  Due  to  noise  from 
vibration  in  the  optical  encoder  mounting,  he  was  unable  to  show  good  agreement  for 
speed  fluctuations  between  measured  results  and  the  model. 

To  the  best  of  the  author’s  knowledge,  no  research  has  been  reported  using  a 
high-fidelity  torsional  model  to  determine  explicit  cylinder-specific  torques  throughout  a 
representative  cycle.  An  engine  crankshaft  displays  torsion  that  varies  during  a  cycle, 
and  along  its  length.  This  twist  absorbs  rotational  energy  that  is  later  released  when  the 
twist  relaxes.  Previous  models  that  consider  the  crankshaft  as  one  rotating  element 
neglect  the  effect  this  twist  produces  on  the  torques  of  individual  cylinders.  But  when  the 
intent  is  to  localize  engine  faults,  it  is  imperative  that  these  differences  between  cylinders 
are  considered. 

C.  OBJECTIVES 

The  primary  objective  of  this  study  is  to  develop  a  method  of  determining 
individual  cylinder  gas  torques  from  measured  time-resolved  angular  positions  at  the  two 
endpoints  of  the  crankshaft,  using  a  high-fidelity  torsional  model  of  the  crankshaft. 
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Calibrated  parameters  will  be  determined  for  the  torsional  model  from  calculated  values 
and  experimental  data. 

Another  objective  is  to  develop  numerical  solution  techniques  for  solving  the 
equations  of  motion  in  both  directions.  Specifically,  a  method  for  direct  integration  of 
the  differential  equations  of  motion  will  be  formulated  that  sets  cyclic  boundary 
conditions  in  the  time  domain.  The  reverse  method,  to  determine  cylinder  gas  torques 
from  tine-resolved  angular  position  data,  must  use  data  from  only  two  measurement 
points  instead  of  all  degrees  of  freedom. 

D.  ORGANIZATION 

Chapter  II  describes  the  development  and  calibration  of  the  crankshaft  torsional 
model.  The  equations  of  motion  for  the  crankshaft  are  derived  and  presented. 

Chapter  III  describes  the  experimental  apparatus  used  to  collect  data  for  this 
study.  Specifications  for  the  test  engine  and  diagrams  for  the  instrumentation  are 
presented.  A  summary  of  the  various  engine  operating  conditions  used  for  the  study  is 
included. 

In  Chapter  IV,  an  explanation  of  the  numerical  methods  used  to  solve  the 
equations  of  motion  is  described.  This  will  include  numerical  methods  for  model 
calibration  as  well  as  cylinder  gas  torque  prediction. 

Chapter  V  shows  the  results  of  the  engine  test  runs  and  data  analysis  that  will 
support  the  thesis  concept. 

Chapter  VI  summarizes  the  study,  presents  conclusions,  and  lists 
recommendations  for  further  research  in  this  area. 
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II.  CRANKSHAFT  TORSIONAL  MODEL 


A.  PHYSICAL  SYSTEM 

The  engine  used  in  this  research  is  a  Detroit  Diesel  (3-53  Series)  three-cylinder 
two-stroke  Diesel.  The  crankshaft  (Figure  1)  is  supported  by  four  main  journal  bearings, 
and  includes  counterweight  lobes  on  four  of  the  six  crankwebs.  The  front  of  the  engine  is 
to  the  left  in  the  diagram.  A  press  fitted  gear  at  the  crankshaft  nose  drives  the  oil  pump, 
and  auxiliary  loads  are  driven  by  the  timing  gear,  located  just  forward  of  the  flywheel 
(not  shown). 


COUNTERWEIGHT 


connecting  rod 
JOURNAL 


REAR  MAIN 
SEARING  JOURNAL 
/ 

/ 


Figure  1.  Crankshaft.  From  Ref  [19]. 


An  idealized  mass-elastic  torsional  model  is  used  to  mathematically  describe  the 

angular  motion  of  the  crankshaft  (Figure  2).  This  model,  originally  developed  by  Hudson 

[Ref  18],  has  been  refined  for  the  present  study.  Specifically,  additional  load  torques 

were  added  to  the  model  to  account  for  the  effect  of  the  oil  pump  and  the  auxiliary  loads, 

constant  parasitic  force  was  used  to  model  the  piston  ring  friction,  and  the  effect  of 

reciprocating  torques  was  added.  The  model  consists  of  six  concentrated  masses 
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connected  by  massless  idealized  shafting.  The  mass  concentrations  are  centered  at  the 
optical  encoder  mounting,  the  three  crankpins,  the  flywheel,  and  the  dynamometer  rotor. 
Torsional  rigidity  and  damping  are  indicated  by  K  and  C,  respectively.  Gas  torques  and 
load  torques  are  applied  at  the  mass  concentrations,  as  indicated  in  parenthesis. 


Cyl#1  Cyl  #2  Cyl#3 


B.  EQUATIONS  OF  MOTION 

Using  the  model  from  Figure  2,  a  set  of  six  second  order  differential  equations  are 
developed  to  describe  the  rotational  dynamics  of  the  crankshaft.  Since  no  part  of  the 
crankshaft  is  fixed,  the  model  requires  six  separate  angular  position  indications.  These 
are  the  crank  angles  at  each  of  the  lumped  mass  points  designated  in  the  model,  and  they 
are  designated  0i  through  06. 
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J\6\  +  C\i{6\  —  0i)  +  K\{6\  —  62)  —  —Tpmp  (l) 

(72 + Jlrejd2 + cm -  ft)  • + K0. 2  -  ft)  ■ +  C2<ft  -  ft)  ■ +  Kli&i  -  ft)  • +  Oft =7^7^ *)-£,,  (2) 

(73  +7*Jft  • +Q<ft-ft)  + J&(ft  -ft)  ■ +C34(ft-ft)  +  ^3(ft-ft)+C3ft  =4>,(0+4«(0  -  '  <3> 

(74 + 74w)& + C^4  -  ft) + K0*  -  ft) + Gs(04  -  ft) + fo(04  -  ft) + Gft =7^,(0 + rVcc(r)  -  Tpar  (4) 

7505  +  C«(ft  -  04)  +  7:4(05  -  04)  +  C 56(05  -  06 )  +  Ks(05  -  06)  =  -F**  (5) 

7  606  +  Cs 6(06  —  05)  +  7l5(06  —  05)  =  —Tload  (6) 

For  simplicity,  the  six  equations  can  be  combined  into  one  matrix  equation: 


or 

[j]e+[C]d+[K]e=fr}  (7) 
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c. 


CALCULATION  OF  MODEL  PARAMETERS 


The  moments  of  inertia  and  the  torsional  rigidities  can  be  analytically  calculated 
for  the  model,  as  described  in  detail  in  Appendices  A  and  B.  It  is  assumed  that  the  values 
for  damping  will  be  very  low.  Final  values  for  the  model  are  listed  in  Table  1 . 


Table  1.  Model  Parameter  Values 


(lbf*in*sec2/rad) 

(106  lbf*in/rad) 

(lbf*in*sec/rad) 

J. 

0.02443 

K, 

3.22 

C\2 

6.01 

h 

0.2482 

k2 

7.00 

C23 

0.01 

h 

0.1462 

k3 

7.00 

C34 

0.01 

h 

0.2482 

K4 

C45 

0.01 

h 

7.2220 

k5 

1.304 

c56 

0.01 

h 

c2 

Jrec 

0.04955(1 -cos29) 

c3 

Jrec,av2 

0.02478 

c4 

Values  for  the  friction  and  auxiliary  loads  are  more  difficult  to  determine 


analytically.  Additionally,  they  will  vary  depending  on  the  load  and  speed  of  the  engine. 
Appendix  A  contains  a  theoretical  analysis  of  friction  and  load  torques,  and  this  analysis 
was  used  to  formulate  an  estimate  of  the  expected  magnitudes  of  friction  and  load 
torques.  The  values  actually  used  in  the  model  are  listed  in  Table  2. 


Table  2.  Model  Friction  and  Load  Values 


RPM 

Load 

(ft*lbf) 

Tload 

(in*lbf) 

Tpmp 

(in*lbf) 

Taux 

(in*lbf) 

Fpardbf) 

(per  cylinder) 

KuZil 

80 

9.5 

160 

98 

9.5 

168 

100 

9.5 

160 

95 

9.5 

300 

155 

9.5 

300 

165 

2000 

160 

1920 

9.5 

460 

210 

Values  for  the  nonlinear  model  parameters  (Trec,  Jrec,  and  Tcy0  are  determined  as 
functions  of  0  or  t.  Their  derivation  is  described  in  detail  in  Appendices  A  and  B. 
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III.  EXPERIMENTAL  METHODS 


A.  ENGINE  DESCRIPTION 

The  Engine  used  in  this  research  was  a  Detroit  Diesel  3-53  Series  engine,  with 
characteristics  listed  in  Table  3.  For  this  study  the  front  of  the  engine  is  designated  as  the 
end  where  the  pulley  would  be  located;  the  rear  is  the  flywheel  end.  The  cylinders  are 
numbered  consecutively  from  front  to  rear,  so  that  cylinder  #1  is  the  farthest  from  the 
flywheel.  Cylinders  are  naturally  aspirated;  a  roots  blower  provides  a  positive  crankcase 
pressure  that  is  proportional  to  engine  speed  [Ref  18].  The  engine  is  considered  to  be  a 
typical  example  of  a  Diesel  engine. 


Table  3.  Engine  Characteristics.  From  Ref  [19] 


Model 

5033-500  IN 

Type 

In-line  two-stroke  compression  ignition 

Number  of  Cylinders 

3 

Number  of  Main  Bearings 

4 

Firing  Order 

1-3-2;  Clockwise  Rotation 

Exhaust  Valves  per  Cylinder 

4 

Displacement  per  Cylinder 

53  ir?  ~1 

Compression  Ratio 

21.0:1 

Bore 

3.875  in 

Stroke 

4.50  in 

Max  Rotation  Speed 

2800  RPM 

Peak  Torque 

198  ft*lbf  @  1500  RPM 

Max  Power  Output 

92BHP 

The  engine  has  been  slightly  modified.  The  front-end  pulley  was  removed  for 
mounting  of  the  optical  encoder  to  the  crankshaft,  and  the  alternator  was  removed  (Figure 
3).  The  engine  was  mounted  on  a  Superflow  engine  test  stand,  and  was  loaded  by  an  SF- 
901  Water  Brake  Dynamometer  (Figure  4). 
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Figure  4.  Engine  Test  Stand  (Side  View) 
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B.  ENGINE  CYCLE  ANALYZER  (ECA) 


The  Superflow  Engine  Cycle  Analyzer  (Figure  5)  is  a  PC  based  data  acquisition 
system.  A  sensor  interface  collects  an  engine  load  signal  from  the  dynamometer,  a  crank 
angle  and  TDC  signal  from  the  optical  encoder,  and  pressure  signals  from  the 
piezoelectric  pressure  transducers  mounted  in  the  glow  plug  sockets  for  each  cylinder. 
This  information  is  passed  to  a  data  acquisition  computer,  which  is  used  to  store  and 
display  the  pressure  data.  Raw  pressure  data  are  collected  and  phase-lock  ensemble 
averaged  over  1 1  cycles,  then  used  in  the  numerical  analysis  programs.  [Ref  20] 


Figure  5.  Instrumentation  Schematic 
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C.  OPTICAL  ENCODER 

A  Heidenhain  incremental  Rotary  Encoder  [Ref  21]  was  used  to  collect  time- 
resolved  angular  position  at  the  crankshaft  nose.  The  encoder  consists  of  a  flat  optical 
disk,  rigidly  attached  to  the  rotating  shaft,  with  a  specified  number  of  evenly  spaced 
windows  etched  near  the  perimeter  (Figure  6).  A  signal  is  generated  by  photoelectric 
scanning  of  the  disk  as  it  rotates.  The  output  signal  is  a  TTL  square  wave,  where  highs 
correspond  to  a  window  passing  in  front  of  the  detector.  Measurement  of  the  leading 
edge  of  the  square  wave  corresponds  to  a  time  stamp  for  a  specific  angular  position.  The 
time  differences  inversely  correspond  to  the  average  speed  of  the  shaft  as  it  rotates 
through  the  incremental  angle.  Encoders  with  720  and  3,600  windows  were  available, 
but  in  either  case  720  counts  per  revolution  were  collected,  for  an  angular  resolution  of 
0.5°.  During  a  run  data  were  collected  for  1 1  cycles  at  the  encoder  for  a  total  of  7,920 
time  stamps. 


Figure  6.  Optical  Disk  Representation 


The  Optical  Encoder  shaft  was  mounted  to  the  end  of  the  crankshaft  with  a 
flexible  coupling  (Figure  7).  The  coupling  allows  for  radial  and  axial  vibration  of  the 
crankshaft  that  would  damage  the  encoder,  because  the  endplay  of  the  crankshaft  exceeds 
the  design  specifications  for  the  optical  encoder  without  the  protective  coupling  installed. 
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Hudson  [Ref  18]  collected  his  data  with  the  encoder  mounted  directly  to  the  crankshaft, 
resulting  in  extremely  noisy  data  that  did  not  compare  favorably  with  predicted  data. 
Additionally,  excessive  crankshaft  radial  vibrations  damaged  several  encoders.  The 
coupling  transmits  the  angular  position  of  the  crankshaft  nose  to  within  an  accuracy  of 
10"  (4.85e-05  radians)  [Ref  21].  An  additional  effect  of  the  coupling  is  a  high  frequency 
torsional  oscillation  due  to  the  natural  frequency  of  the  coupling/rotor  combination.  This 
is  discussed  in  detail  in  Appendix  C,  and  was  not  a  significant  problem  since  the  signal  of 
interest  was  at  a  much  lower  frequency. 


Figure  7.  Optical  Encoder  and  coupling.  From  Ref  [21] 

The  body  of  the  encoder  is  rigidly  mounted  to  the  engine  block  to  minimize  noise 
due  to  vibration  (Figure  8).  The  mounting  of  the  encoder  is  extremely  important. 
Hudson  used  several  different  mounting  schemes,  before  settling  on  the  mount  used  again 
in  this  study,  which  works  very  well  to  ensure  engine  vibration  does  not  affect  the 
encoder  measurements. 
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Figure  8.  Optical  Encoder  Mounting.  From  Ref  [18] 


D.  MAGNETIC  INDUCTION  PROXIMETER 

A  Bentley  Nevada  3000  series  190  Proximitor  system  was  used  to  detect  passage 
of  the  flywheel  ring  gear  teeth.  The  system  consists  of  a  ferromagnetic  eddy  current 
detector,  which  outputs  a  negative  voltage  that  is  a  function  of  the  distance  between  the 
probe  end  and  a  ferromagnetic  surface.  A  TTL  conversion  circuit  triggers  a  step  change 
in  voltage  when  the  proximeter  output  exceeds  a  certain  level,  corresponding  to  a 
distance  of  about  XA  inch.  The  probe  was  mounted  on  a  bracket  fastened  across  the  edge 
of  the  flywheel  access  panel,  so  that  it  saw  the  sides  of  the  gear  teeth  as  they  passed 
(Figure  9).  The  output  of  the  circuit  was  a  square  wave  TTL  signal;  the  leading  edge  of 
each  wave  corresponding  to  the  passage  of  a  gear  tooth  beneath  the  probe.  The  TTL 
output  signal  was  sent  directly  to  the  MDA  for  data  collection.  There  are  126  teeth  on  the 
flywheel  ring  gear,  and  during  a  run  data  were  collected  for  63  cycles,  for  a  total  of  7938 
time  stamps. 
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Figure  9.  Proximeter  mounting 


E.  MODULATION  DOMAIN  ANALYZER 

A  Hewlett  Packard  533 10A  Modulation  Domain  Analyzer  (MDA)  was  used  to 
collect  the  time  stamp  data  from  the  optical  encoder  and  the  proximeter.  In  either  case,  a 
TDC  indicator  (an  output  from  the  ECA)  triggered  the  MDA.  It  received  a  TTL  signal 
and  recorded  a  time  stamp  at  the  leading  edge  of  each  wave.  The  MDA  was  able  to 
collect  up  to  8,000  data  points  at  a  time.  A  single  ran  collected  1 1  cycles  from  the  optical 
encoder  or  63  cycles  from  the  flywheel  proximeter,  as  described  previously.  A  data 
acquisition  computer  controlled  the  MDA  and  received  a  transferred  file  containing  the 
time  stamps. 


F. 


TEST  MATRIX 


A  series  of  data  runs  were  performed  to  test  the  validity  and  consistency  of  the 
model  at  varying  engine  speeds  and  loads.  During  a  run  data  were  collected  for  rotational 
speed  at  the  flywheel  and  crankshaft  nose,  cylinder  pressures  for  the  three  cylinders, 
dynamometer  load,  and  atmospheric  pressure.  A  data  series  was  composed  of  data 
collected  during  a  single  run  of  the  engine,  at  varying  loads  and  speeds,  normally  taking 
about  an  hour  on  a  single  day.  A  prefix  letter,  such  as  “S,”  designated  a  particular  series 
so  that  comparisons  could  be  restricted  to  data  taken  during  a  single  operation  of  the 
engine.  This  was  intended  to  eliminate  any  variations  in  operation  that  might  take  place 
as  the  engine  condition  varied  over  time.  Table  4  shows  the  elements  of  each  data  series, 
indicating  what  variations  in  load  and  speed  make  up  each  one.  Table  5  lists  specific 
information  for  each  data  run. 


Table  4.  Speeds  and  Loads  for  Data  Series 


ENGINE 

DYNAMOMETER  LOAD  (FT*LBF) 

SPEED  (RPM) 

80 

100 

135 

160 

180 

1000 

s,u,v,w 

S,T,U,V,W 

S,T,U,V,W 

1500 

S,T,U,V,W 

S,T,U,V,W 

2000 

T,U,V,W 

V 

Table  5.  Data  Run  Information 


DATA  RUN 

Date 

Ambient 

Pressure 

(psia) 

Comments 

S 

02  Sep  98 

14.593 

Heli  -  Cal  coupling  used;  slippage  of 
encoder  shaft 

T 

1 1  Sep  98 

14.662 

Heli  -  Cal  coupling  used;  TDC  lag 

U 

08  Oct  98 

14.819 

Heidenhain  K17  coupling  used;  TDC  lag 

V 

14  Oct  98 

14.730 

Heidenhain  K17  coupling  used;  TDC  lag 

W 

21  Oct  98 

14.706 

Heidenhain  K17  coupling  used;  no  lag 
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IV.  METHODS  FOR  CYLINDER  PRESSURE  AND  TORQUE  PREDICTION 

The  torsional  model  may  be  used  to  predict  the  motions  of  the  crankshaft  from 
measured  applied  forces,  or  it  may  be  used  to  predict  the  forces  from  the  measured 
motions.  In  order  to  do  this,  the  differential  equations  of  motion  must  be  solved  in  both 
directions.  First,  two  methods  will  be  described  for  determining  the  angular  positions  0i 
through  06,  given  the  measured  gas  torques  from  the  three  cylinders.  These  solution 
methods,  referred  to  as  direct  integration  methods,  will  be  used  to  test  the  validity  of  the 
torsional  model  and  calibrate  the  parameters.  Second,  a  method  will  be  described  for 
determining  the  individual  cylinder  gas  torques  T|cyl  through  T3Cyi,  given  the  time- 
resolved  angular  position  at  the  two  ends  of  the  crankshaft,  0i  and  @5.  This  final  solution 
method,  called  the  inverse  method,  will  be  used  for  detection  of  cylinder  faults  from 
measurements  of  crankshaft  rotational  velocity. 

A.  PREDICTION  OF  PHASE  DEVIATION  FROM  CYLINDER  TORQUES 

The  equations  of  motion  for  the  crankshaft  model  constitute  a  system  of  non¬ 
linear,  second-order  ordinary  differential  equations  (ODEs).  Calibration  of  the  values  for 
the  torsional  model  will  be  conducted  by  first  solving  the  differential  equation  for  {0}, 
given  the  cylinder  indicated  torques.  The  predicted  phase  deviation  and  twist  determined 
from  {0}  can  be  compared  to  measured  values  to  determine  validity  of  the  model. 

1.  Time-marching  O.D.E.  method 

The  first  method  is  a  direct  integration  of  the  ODEs  in  the  time  domain.  This  is 
accomplished  numerically  using  a  fourth-  and  fifth-order  Runge-Kutte  method.  First,  the 
six  second-order  ODEs  must  be  converted  to  12  first-order  ODEs  as  follows: 
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di=0\ 

4=4 

6q  =4 

^lO  =  4 

0i  i = 4 

012  =  4 

■A 4  —  C|2  08  +CJ2^7  — -^1^2  ■*'•^101  ~~Tpmp 
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(8) 

These  equations  cannot  be  solved  all  in  one  step,  however,  because  they  are 
nonlinear  due  to  the  dependence  of  Jrec  on  0,  and  because  the  TcyiS  and  Trecs  are  functions 
of  time.  Instead,  a  “time  marching”  method  is  used  where  the  equations  are  solved  over 
small  steps  and  the  final  condition  of  each  step  becomes  the  initial  condition  for  the 
subsequent  step.  The  values  of  Tcyi,  Trec,  and  Jrec  can  then  be  approximated  as  a  constant 
value  over  the  step,  or  as  a  linear  interpolation  within  the  function  describing  the 
equation. 

This  method  requires  significant  computing  time.  This  is  because  a  series  of 
“shooting”  iterations  must  be  conducted  to  determine  the  initial  conditions  that  yield 
cyclic  conditions  for  the  representative  cycle  (i.e.,  a  periodic  solution).  Since  the  domain 
is  an  assumed  representative  cycle,  it  follows  that  the  values  of  6  and  6  at  the  end  of  the 
cycle  must  match  those  at  the  beginning  of  the  cycle.  The  initial  angular  velocities 
chosen  will  have  a  significant  effect  on  whether  or  not  the  solution  0  is  “cyclic.” 


20 


2. 


Finite  Element  method 


Although  the  integration  of  the  differential  equations  is  an  initial  value  problem, 
because  it  is  assumed  to  describe  a  cyclic  process  the  solution  of  the  equation  is  known  at 
future  times.  That  is,  at  the  end  of  one  representative  cycle,  we  expect  all  the  angular 
positions  to  be  increased  by  exactly  2%  radians.  An  alternative  method  for  solution 
avoids  the  shooting  iterations  by  setting  “boundary”  conditions  in  time,  instead  of  initial 
conditions.  The  problem  is  then  treated  as  a  boundary  value  problem  in  time,  and  a  Finite 
Element  Method  (FEM)  is  developed  to  solve  a  second-order  differential  equation  for  0 
as  a  function  of  t.  This  method  is  based  on  Kwon  and  Bang  [Ref  22], 

The  weak  formulation  of  the  weighted  residual  method  is  used  to  approximate  the 
solution  to  a  second-order  matrix  differential  equation.  To  accomplish  this,  the  weighted 
average  of  the  residual  over  the  domain  is  set  to  zero: 

I  =  J  W  +  [CW  +  [K]6  -  {r}}?r  =  {o}  (9) 

'() 

and  then  simplified  to: 

'/  '/  lf 

\{w[c]-w[j}e  +  w[Klp}lt=jwfr}dt-\w[jlp}  (10) 

*n  r0 

where  w  is  the  weighting  function,  and  0  is  a  vector  corresponding  to  the  angular  position 
at  the  six  degrees-of-freedom. 

.  A  Galerkin  Finite  Element  formulation  is  developed,  using  the  sum  of  simple 
piecewise  linear  shape  functions  to  approximate  the  more  complex  real  function.  The 
shape  functions  used  are  set  up  to  be  1  at  a  node,  linearly  decreasing  to  0  at  the  adjacent 
nodes.  The  value  of  the  function  at  any  point  is  approximated  as  a  linear  combination  of 
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the  values  of  the  two  shape  functions  for  the  adjacent  nodes.  Figure  10  shows  these 
linear  shape  functions  for  a  hypothetical  element.  The  function  T(t)  is  approximated  as 
T(t)  =  Hi(t)T(t|)  +  H2(t)T(tr)  between  the  nodes.  This  becomes  the  trial  function  for 
Galerkin’s  method,  and  the  test  functions  are  W]=Hi(t)  and  W]=H2(t). 


T(tl) 


4r~  h - ^ 


HIT 


Figure  10.  Linear  Shape  Functions.  After  Ref  [22] 


For  a  one-dimensional  function,  the  shape  functions  are  determined  by: 

H\  (0  =  -~7~  h  =  tR-tL  •  (11) 

h  h 

But  when  used  for  a  6x6  matrix  equation,  the  test  functions  must  also  be  expanded  into 
the  following  matrices: 

[tf,  ]  =  ,[/] . [H2]  =  H2[I] . [77,  ]  =  -  j[7] . [H2]  =  |[7]  (12) 

h  h 

When  substituted  in  Equation  (10)  and  solved  over  the  element  from  tL  to  tR,  the  result  is: 


where 


) 


(14) 


and 
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Here,  the  element  matrix  [Ke]  is  12x12  and  the  element  torque  vector  { F0}  is  12x1.  The 
values  of  the  torque  at  the  two  nodes,  {TL}  and  {TR},  are  6x1  vectors.  Equation  (13)  may 
be  solved  for  the  6x1  vectors  {0L}  and  {0R},  which  are  the  approximate  solutions  at  the 
nodes.  The  solution  used  in  this  study  has  720  elements  and  721  nodes  to  solve  for  one 
representative  cycle.  The  element  matrices  for  each  node  are  assembled  into  a 
4326x4326  finite  element  matrix  [Kfe]  and  the  element  torque  vectors  are  assembled  into 
a  4326x1  finite  element  torque  vector  {Ffe}.  Boundary  conditions  are  established  by 
defining  the  value  of  the  6x1  vector  {0}  at  the  end  nodes.  This  is  accomplished  by 
setting  the  first  and  last  six  lines  of  the  finite  element  matrix  [Kfe]  to  identity,  and  setting 
the  first  six  and  last  six  values  in  the  finite  element  torque  vector  {Ffe}  to  the  boundary 
values.  The  solution  {0}  (a  4326x1  vector,  the  720  6x1  nodal  solutions  {0;}  stacked 
vertically)  is  then  found  from  the  matrix  equation: 

[Kt'K>'s={F“}  <i6) 

This  method  shows  a  marked  improvement  in  computing  efficiency  over  the  time¬ 
marching  method.  In  addition  to  being  about  three  times  as  fast  for  each  program  ran,  the 
Finite  Element  Method  avoids  the  shooting  iterations  required  for  the  time-marching 
method,  which  had  to  be  repeated  as  many  as  five  times  for  each  solution.  A  comparison 
of  the  Phase  Deviation  found  with  both  methods  is  shown  in  Figure  11.  The  Phase 
Deviation  Plot  will  be  described  in  detail  in  the  next  chapter. 


B.  PREDICTION  OF  TORQUES  FROM  MEASURED  SHAFT  SPEEDS 


1.  Solution  of  Matrix  Equation 
a.  Two-Stroke  Engines 

Given  a  properly  calibrated  model,  the  time  dependent  values  of  gas 
torque  for  the  three  cylinders  (Ticyl,  T2cyl,  and  T3cyl)  can  be  determined  from  the  six 
second-order  ODEs  (Equations  1-6).  Figure  12  outlines  the  solution  method.  By 
measurements  at  the  flywheel  and  crankshaft  nose,  the  values  for  0i  and  65  can  be 
determined  for  a  representative  cycle.  Numerical  differentiation  of  this  data  yields  the 
velocities  and  accelerations  at  these  two  points. 

With  properly  calibrated  parameters  for  the  torsional  model,  Equation  (6) 
can  be  solved  for  06,  then  Equation  (5)  can  be  solved  for  04,  and  Equation  (1)  can  be 
solved  for  02.  Now  there  are  four  unknowns  left  (03,  Tlcy|,  T2cyi,  and  T3cyi)  and  three 
equations  (Equations  2,  3,  and  4).  However,  because  this  is  a  two-stroke  engine,  each 
cylinder  is  at  reference  pressure  for  one  third  of  each  cycle,  while  the  exhaust  and  intake 
ports  are  uncovered.  Therefore,  at  any  one  time  during  a  representative  cycle,  one  of  the 
three  cylinder  torques  is  known,  leaving  three  unknowns  and  three  equations.  For 
example,  for  the  first  120  degrees  of  the  cycle,  cylinder  #2  ports  are  open,  so  the  pressure 
in  cylinder  #2  is  reference  pressure  (See  Figure  13).  For  the  first  120  degrees,  T2cyi  can 
be  calculated  from  this  reference  pressure.  Then  03  can  be  calculated  from  Equation  (3), 
and  Ticyi  and  T3cyi  can  be  determined  from  Equations  (2)  and  (4).  The  other  two-thirds  of 
the  representative  cycle  are  solved  in  a  similar  manner. 
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90 


Crank  Anale  (decrees’) 


Figure  13.  Cylinder  Gas  Pressures  (1000  RPM,  100  Ft*lbf) 

Although  the  ports  are  not  actually  uncovered  for  the  first  and  last  10  degrees, 
measurements  show  the  assumption  of  reference  pressure  is  reasonable  for  the  entire  120 
degrees.  The  assumptions  used  here  limit  the  feasibility  of  this  solution  method  to  two- 
stroke  engines  with  three  or  fewer  cylinders. 

b.  Four-Stroke  Engines 

For  a  four  stroke  engine,  two  full  rotations  of  the  crankshaft  must  be 
considered  to  cover  the  power,  exhaust,  intake,  and  compression  strokes.  Over  the 
representative  two  rotation  (one  cycle)  period,  a  cylinder’s  pressure  can  be  assumed  to  be 
equal  to  intake  manifold  pressure  during  the  intake  stroke,  and  exhaust  back  pressure 
during  exhaust  stroke.  Therefore,  for  a  particular  cylinder,  torque  is  known  for  half  of  the 
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cycle.  The  solution  method  above  is  feasible  for  four-stroke  engines  with  four  or  fewer 
cylinders.  A  four  cylinder  engine  would  require  a  seven  degree  of  freedom  model,  which 
could  be  explicitly  solved  as  above,  using  the  assumptions  for  intake  and  exhaust  stroke 
pressure. 

c.  Additional  Assumptions  for  Multiple  Cylinder  Engines 

Engines  with  more  than  four  cylinders  can  be  analyzed  by  this  method  if 
further  assumptions  are  made.  For  instance,  the  cylinder  compression  stroke  can  be 
estimated  as  a  polytropic  compression  of  an  ideal  gas.  For  a  large  engine  with  cylinders 
that  fire  simultaneously,  the  two  cylinders  could  be  lumped  and  considered  as  one  inertial 
mass.  Also,  a  measurement  device  may  be  placed  internal  to  the  engine  to  measure 
angular  velocity  at  a  third  degree-of-freedom. 

2.  Interpolation  of  Data 

The  rotational  speed  data  that  are  collected  at  the  optical  encoder  and  the  flywheel 
consists  of  information  that  is  uniformly  spaced  in  the  angular  position  domain,  not  in  the 
time  domain.  This  arises  because  of  the  nature  of  the  data  collection  (See  Section  III.C 
and  III.D).  The  raw  data  for  0i(t)  and  0s(t)  are  converted  to  values  which  are  evenly 
spaced  in  the  time  domain  for  further  numerical  analysis  (specifically,  this  is  useful  for 
filtering;  see  section  IV.B.3).  This  requires  interpolation  of  the  raw  data.  Interpolation  is 
accomplished  numerically  using  a  cubic  spline  method.  This  means  that  the  curve  is 
assumed  to  be  a  3rd  order  polynomial  with  continuous  slope  at  each  of  the  data  points. 
Interpolated  values  of  0(t)  are  determined  for  a  selected  evenly-spaced  time  basis. 
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3.  Signal  Filtering  by  Fast  Fourier  Transform 

The  measured  data  for  position,  0(t),  and  angular  velocity,  ^-1 1  contain  high 

at 

frequency  components  (where  “high  frequency”  refers  to  frequencies  much  higher  than 
about  three  times  the  rotational  speed).  These  frequency  components  are  due  to  high 
frequency  torsional  vibration  of  the  crankshaft  at  the  various  natural  frequencies,  and 
random  noise  from  unknown  sources.  For  solution  of  the  equations  necessary  to  predict 
torques,  the  position  data  must  be  differentiated  once  to  determine  angular  velocity  and 
twice  to  determine  angular  acceleration  (see  section  IV.B.l).  If  raw  data  were  used  in  the 
analysis,  the  high  frequency  components  would  be  greatly  amplified  by  subsequent 
differentiation.  However,  the  torques  of  primary  interest  in  this  problem  oscillate  at 
about  three  times  the  rotational  velocity  of  the  engine.  Therefore,  the  much  higher 
frequency  components  are  filtered  out  before  differentiation  in  order  to  increase  the 
signal-to-noise  ratio  in  solving  the  equations  for  torque. 

A  phase  deviation  signal  is  derived  from  the  raw  data  by  comparing  it  to  the  mean 
rotational  speed.  This  phase  deviation  signal  is  then  filtered  numerically  using  a  Fast 
Fourier  Transform  (FFT),  removing  the  high  frequency  data,  and  then  performing  an 
inverse  FFT  to  achieve  the  desired  filtered  phase  deviation  signal.  The  cut-off  frequency 
typically  used  was  between  six  and  nine  times  the  rotational  speed  of  the  engine. 

4.  Numerical  Differentiation 

Once  the  angular  position  data  has  been  interpolated  and  filtered,  it  must  be 
differentiated  once  to  obtain  angular  velocity  and  twice  to  obtain  angular  acceleration  as 
functions  of  time.  A  central  difference  technique  is  used,  where  the  numerical  derivative 

at  a  point  is  the  sum  of  the  two  adjacent  differences  divided  by  twice  the  time  difference. 
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The  endpoints  are  exceptions;  they  are  determined  by  the  single  adjacent  difference 
divided  by  the  time  difference  (essentially  a  forward  difference  for  the  first  point  and  a 
backward  difference  for  the  last  point).  No  further  filtering  is  required  after 
differentiation. 

C.  TEST  OF  SOLUTION  METHODS 

The  methods  previously  discussed  are  used  to  solve  the  equations  of  motion  in 
both  directions.  Using  a  set  of  measured  cylinder  indicated  pressures,  the  accuracy  of  the 
numerical  methods  can  be  tested.  The  time-marching  ODE  method  was  used  to  solve  for 
the  crankshaft  time-resolved  angular  positions  from  the  measured  cylinder  indicated 
torques.  The  predicted  0i  and  65  values  were  then  used  in  the  inverse  solution  method  to 
predict  the  cylinder  indicated  torques.  Comparison  of  the  resulting  predicted  torques  to 
the  original  measured  torques  can  be  used  to  quantify  the  accuracy  of  the  numerical 
solution  methods.  The  test  results  for  individual  cylinder  gas  torques  are  shown  in  Figure 
14,  and  the  results  for  total  gas  torque  is  shown  in  Figure  15.  These  results  show  that  the 
numerical  methods  tend  to  introduce  a  2%  peak-to-peak  error  and  a  4  degree  lag  in  the 
predicted  total  gas  torque. 
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Test  Run 


Figure  15.  Test  of  Numerical  Methods  for  Total  Torque 
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V.  RESULTS 


A.  CALIBRATION  OF  INERTIAL  MODEL 

Although  the  stiffness  and  inertial  parameters  of  the  equations  of  motion  can  be 
analytically  estimated,  some  fine-tuning  of  the  values  is  required  to  ensure  the  model 
matches  the  actual  crankshaft.  Given  the  measured  torques  in  the  three  cylinders,  and  the 
methods  of  the  previous  chapters,  the  equations  of  motion  may  be  solved  for  0i  through 
06.  Experimental  data  is  collected  for  the  two  measurement  points,  0i  and  05.  Two 
useful  comparisons  between  the  measured  and  predicted  time-resolved  angular  positions 
are  the  Phase  Deviation  Plot  and  the  Crankshaft  Twist  Plot. 

The  Phase  Deviation  Plot  shows  the  oscillation  of  the  angular  position  about  a 

theoretical  mean  rotating  position.  It  is  calculated  as  follows: 

£(t)  =  6(t)-mt  (17) 

where  a  is  the  mean  rotational  velocity.  This  is  the  same  comparison  plot  used  by 

Hudson  [Ref  18].  For  a  shaft  rotating  at  a  steady  angular  velocity,  the  phase  deviation 

would  be  zero.  The  measured  phase  deviation  shows  the  crankshaft  position  advancing 

during  the  power  stroke  of  each  cylinder,  then  retreating  during  the  subsequent 

compression  of  the  next  cylinder.  A  comparison  of  the  measured  and  predicted  phase 

deviation  is  shown  for  the  crankshaft  nose  and  the  flywheel  (Figure  16).  Plots  for  other 

runs  are  found  in  Appendix  D.  The  phase  deviation  plot  is  particularly  sensitive  to 

assumed  values  of  friction  and  load  torque,  and  is  also  useful  for  validating  the  model 

inertias.  For  the  1000  RPM  100  Ft*lbf  data,  phase  deviation  shows  a  maximum  error  of 

about  9%  at  the  peaks,  which  correspond  to  the  cylinder  compression  strokes.  Generally, 

the  model  predicted  phase  deviation  is  within  5%  of  the  measured  phase  deviation  for 
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both  the  flywheel  and  the  crankshaft  nose.  The  errors  are  calculated  as  a  percentage  of 
the  maximum  variation  in  phase  deviation. 

The  Crankshaft  Twist  Plot  is  simply  the  value  of  0i-05  as  a  function  of  time.  It  is 
a  measure  of  the  total  twist  along  the  entire  length  of  the  crankshaft.  A  comparison  of 
measured  and  predicted  crankshaft  twist  is  shown  in  Figure  17.  This  plot  is  particularly 
useful  for  validating  the  magnitudes  of  the  torsional  rigidities  used  in  the  model.  The 
peak  values  of  the  twist  occur  at  the  peaks  of  the  gas  torque,  during  the  power  stroke  for 
each  cylinder.  As  expected,  the  amount  of  twist  is  largest  for  cylinder  #1,  the  farthest 
from  the  flywheel,  and  successively  lower  for  cylinders  #2  and  #3.  Correct  values  of 
torsional  rigidity  for  the  model  should  result  in  predicted  twists  comparable  to  the 
measured  value.  Comparison  of  the  model  predicted  twist  and  the  measured  twist  is 
shown  in  Figure  17,  with  additional  plots  in  Appendix  D.  A  max  error  of  about  17%  is 
seen  at  the  peak  twist  values,  corresponding  to  the  cylinder  power  strokes.  Generally,  the 
model  predicted  twist  is  within  about  5%  of  the  measured  twist. 

As  discussed  in  Appendix  C,  analysis  of  the  crankshaft  natural  frequencies  can 
also  be  used  to  validate  the  model  parameters.  The  analytical  values  derived  for  the 
crankshaft  inertias  are  considered  to  be  reasonably  accurate,  so  the  only  parameters  to  be 
adjusted  are  the  torsional  rigidities.  These  are  set  by  the  comparison  of  measured  and 
predicted  natural  frequencies  and  modes  as  discussed  in  Appendix  C.  Further  arbitrary 
adjustment  of  the  parameters  to  correct  the  errors  in  the  Phase  Deviation  and  Crankshaft 
Twist  plots  is  not  supported. 
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B.  COMPARISON  OF  PREDICTED  AND  MEASURED  TORQUES 


Inverse  solution  of  the  equations  of  motion  allows  prediction  of  the  cylinder  gas 
torques,  given  the  time-resolved  angular  position  at  two  points,  0i  and  65.  Comparison  of 
measured  and  predicted  gas  torques  for  the  individual  cylinders  is  shown  in  Figure  18. 
Comparison  is  not  good  for  individual  torques.  Besides  quantitative  errors  of  over  50% 
at  certain  points,  the  predicted  torques  show  misplaced  and  inappropriate  peaks. 
However,  it  appears  that  the  errors  for  a  particular  predicted  cylinder  torque  have 
corresponding  offsetting  errors  in  the  predicted  torque  for  the  other  cylinders.  This  is 
evident  when  measured  and  predicted  total  gas  torques  are  compared  (Figure  19).  The 
predicted  total  gas  torque  shows  peak-to-peak  errors  of  less  than  5%,  plus  a  phase  lag  of 
about  5-15  degrees.  Some  of  this  error  is  from  the  numerical  solution  methods,  as 
discussed  in  the  previous  chapter.  This  agreement  is  good  enough  to  be  used  for 
localized  fault  detection.  For  this  particular  run  of  the  engine,  cylinder  #1  gas  pressure  is 
significantly  lower  than  the  other  two  cylinders,  and  the  predicted  results  detect  this 
anomaly.  Plots  for  the  other  data  runs  are  contained  in  Appendix  D. 
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VI.  SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 


A.  SUMMARY 

A  three  cylinder,  two  stroke  Diesel  engine  was  instrumented  with  a  proximeter 
and  an  optical  encoder  for  time-resolved  angular  position  measurement  at  the  flywheel 
and  crankshaft  nose.  A  torsional  model  for  the  engine  crankshaft  was  developed  and  the 
corresponding  equations  of  motion  were  formulated.  Two  separate  numerical  solution 
methods  were  developed  to  solve  for  the  angular  positions,  given  the  measured  cylinder 
gas  torques.  These  methods  were  used  to  calibrate  the  parameters  of  the  torsional  model. 
An  inverse  solution  method  was  devised  to  determine  the  cylinder  gas  torques,  given  the 
time-resolved  angular  positions  at  two  of  the  degrees  of  freedom;  the  flywheel  and  the 
crankshaft  nose.  This  inverse  solution  method  was  shown  to  be  applicable  for  two-stroke 
engines  of  three  or  fewer  cylinders,  or  for  four-stroke  engines  of  four  or  fewer  cylinders. 
The  predicted  cylinder  gas  torques  were  compared  to  measured  cylinder  gas  torques. 

B.  CONCLUSIONS 

The  torsional  model  accurately  describes  the  dynamics  of  the  actual  crankshaft. 
Experimental  data  demonstrated  that  the  model  correctly  predicted  phase  deviation  at  the 
crankshaft  endpoints  with  an  error  of  less  than  5%.  The  model  predicted  crankshaft  twist 
with  an  error  of  less  than  20%.  Predicted  natural  frequencies  from  the  model  agreed  with 
the  measured  frequency  spectrum  to  within  5  Hz  for  the  three  vibration  modes  observed. 

The  Finite  Element  Method  (FEM)  for  direct  integration  of  the  equations  of 
motion  agreed  with  the  Time-marching  ODE  method  to  within  1%,  and  it  reduced 
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computational  time  by  a  factor  of  100  because  no  iterations  were  necessary.  It  was  used 
as  the  primary  means  for  direct  integration  of  the  equations  of  motion  in  this  study. 

The  inverse  method  for  predicting  cylinder  gas  torques  showed  significant  errors 
in  predicted  individual  cylinder  gas  torques.  Quantitative  errors  of  over  50%,  as  well  as 
significant  wave  shape  errors,  make  this  method  inadequate  for  reliable  prediction  of 
individual  cylinder  torques.  However,  it  is  the  author’s  opinion  that  this  error  originates 
in  the  numerical  method  used.  Specifically,  signal  filtering  tends  to  create  errors  in  the 
endpoints  for  the  representative  cycle.  Further  analysis  of  this  problem  may  result  in 
successful  prediction  of  individual  cylinder  torques. 

The  inverse  method  is  successful  in  predicting  total  cylinder  gas  torque. 
Predicted  total  gas  torque  errors  were  less  than  5%,  with  slight  phase  lag  errors  of  5-15 
degrees.  The  predicted  total  gas  torque  successfully  detected  a  low  pressure  in  cylinder 
#1,  showing  that  the  method  is  capable  of  localizing  certain  faults  to  a  particular  cylinder. 

C.  RECOMMENDATIONS 

The  failure  of  the  inverse  method  to  predict  individual  cylinder  torques  is  most 
likely  due  to  problems  with  the  numerical  solution  methods.  The  FFT  signal  filtering 
induces  some  errors,  which  were  not  sufficiently  corrected.  The  engine  speed  was  not 
exactly  steady  during  data  collection,  so  the  filtering  process  backs  out  a  monotonic  trend 
in  the  phase  deviation.  That  is,  the  filtered  phase  deviation  is  no  longer  exactly  cyclic.  A 
correction  of  some  sort  should  be  made  for  this  linear  error.  The  process  of  signal 
filtering  also  tends  to  alter  the  endpoints  of  the  representative  cycle,  inducing  significant 
errors  in  the  calculated  torques  at  the  endpoints.  An  alternative  signal  filtering  process, 
which  avoids  these  errors,  might  correct  the  errors  in  the  results. 
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The  assumption  of  constant  parasitic  force  to  model  piston  ring  friction  is  not 
validated.  The  presence  of  crank-angle  specific  friction  "sticking"  points  would  have  a 
significant  effect  on  the  results.  One  solution  would  be  to  use  a  motoring  dynamometer 
on  the  engine  to  produce  a  graph  of  the  engine  friction  as  a  function  of  crank  angle.  Also, 
a  more  detailed  analysis  of  theoretical  piston  ring  friction  could  lead  to  more  accurate 
modeling. 

From  measurements  obtained  in  this  study,  there  is  an  unknown  fault  causing 
cylinder  #1  to  have  a  lower  gas  pressure  than  the  other  two  cylinders.  As  a  first  step,  the 
fuel  injectors  for  cylinders  #1  and  #2  should  be  swapped  in  order  to  determine  the  cause 
of  the  low  pressure  in  cylinder  #1.  For  follow-on  experimentation,  data  should  be 
collected  for  engine  runs  with  known  faults.  For  instance,  a  defective  fuel  injector  could 
be  installed  to  test  the  method’s  ability  to  detect  a  specific  fault. 

The  use  of  angular  speed  measurement  internal  to  the  engine  would  expand  this 
method  to  engines  with  more  cylinders.  Although  this  would  be  a  difficult  process  for  an 
existing  engine,  mass-produced  engines  might  have  such  an  internal  detector  installed  for 
relatively  little  extra  cost. 

The  method  for  determining  TDC  for  cylinder  #1  is  inadequate.  Currently,  the 
procedure  of  Appendix  F,  Ref  [18]  is  used  to  orient  the  TDC  signal  on  the  encoder  to 
TDC  for  the  engine.  But  TDC  for  the  engine  is  established  by  a  mark  inscribed  on  the 
crankcase  and  the  forward  counterweight  on  the  cam  follower  shaft.  Although  this  shaft 
is  directly  geared  to  the  crankshaft,  gear  backlash  results  in  an  error  of  one  or  two  degrees 
when  the  engine  is  rotated  to  TDC.  This  small  error  has  a  significant  effect  on  the 
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magnitudes  of  the  cylinder  indicated  torques.  A  better  method  would  be  a  mark  inscribed 
on  the  flywheel  and  shroud  to  ensure  correct  TDC  alignment. 

A  sprocket  and  proximeter  assembly  might  be  a  more  useful  means  of  collecting 
data  at  the  front  end  of  the  crankshaft.  A  42  tooth  sprocket  has  been  obtained  which  may 
be  mounted  on  the  crankshaft  nose  with  the  pulley  mounting  bolt.  Since  the  number  of 
teeth  on  the  flywheel  (126)  is  a  whole  number  multiple  of  42,  a  precise  alignment  could 
be  made  to  calibrate  the  static  phase  difference  between  the  two  ends  of  the  crankshaft  to 
zero.  Then  the  instantaneous  twist  of  the  crankshaft  could  be  very  accurately  measured 
during  operation,  similar  to  the  method  described  by  Mauer  and  Watts  [Ref  2]. 
Additionally,  this  data  collection  method  would  eliminate  the  natural  frequency  torsional 
vibrations  of  the  optical  encoder  and  flexible  coupling. 
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APPENDIX  A.  GEOMETRY  OF  ROTATING  COMPONENTS 


While  a  mass-elastic  model  has  already  been  developed  for  this  crankshaft  [Ref 
23],  it  is  necessary  to  independently  calculate  the  model  parameters  in  order  to  verify 
them,  and  to  account  for  differences  in  the  specific  engine  used  in  the  research.  This 
analysis  was  carried  out  using  methods  and  equations  from  Wilson  [Ref  24].  The 
descriptions  and  values  for  certain  constants  used  in  subsequent  equations  are  listed  in 
Table  6. 


Table  6.  Equation  Constants 


Svmbol 

Description 

Value 

Deo 

Diameter,  crankpin 

2.50  in 

Djb 

Diameter,  journal  bearing 

3.00  in 

s 

Gravitational  Acceleration 

386  in/sec2 

Kgyr 

Radius  of  Gyration 

Le 

Length,  effective 

Let) 

Length,  crankpin 

1.60  in 

Lib 

Length,  journal  bearing 

1.50  in 

R 

Crankpin  eccentricity 

2.25  in 

Ri  ^ 

Counterweight  Inner  Radius 

1.80  in 

Ro 

Counterweight  Outer  Radius 

4.02  in 

Tct 

Counterweight  Thickness 

0.83  in 

Twb 

Crankweb  Thickness 

1.00  in 

w 

Weight 

Wwb 

Crankweb  width 

3.84  in 

p 

Specific  weight  of  steel 

0.283  lbf/inJ 

1.  Rotating  Inertia 

An  arbitrary  objects  mass  polar  moment  of  inertia  is  calculated  as: 


g 


where  W  is  the  weight  of  the  object,  Kgyr  is  the  radius  of  gyration,  and  g  is  acceleration 
due  to  gravity. 
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The  polar  moment  of  inertia  for  a  main  journal  bearing  is  determined  as  for  a 


As  seen  in  Figure  20,  the  Crankwebs  are  of  three  distinctive  types.  Crankwebs  of 
type  (a)  connect  crankpins  for  cylinders  1  and  3  to  the  outer  journal  bearings,  and 
crankwebs  of  type  (b)  connect  those  same  crankpins  to  the  inner  journal  bearings. 
Crankwebs  of  type  (c)  connect  the  crankpin  for  cylinder  2  to  the  inner  journal  bearings. 


(a)  (b)  (c) 

Figure  20.  Crankweb  Forms.  From  Ref  [25] 


For  all  three  types,  the  core  geometric  shape  is  an  ellipse  with  one  focus  centered  on  the 

journal  bearing  and  the  other  focus  centered  on  the  crankpin.  From  Table  10.13  of 

Wilson  [Ref  24],  the  radius  of  gyration  for  an  ellipse  can  be  calculated  as: 
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(A4) 


_  (a2  +fr2)  2 

ellipse  ~  |  ^  +C 


where  a  and  b  are  the  minor  and  major  axis,  respectively,  and  c  is  the  offset  of  the  center 
of  gravity  from  the  axis  of  rotation.  J  can  then  be  calculated  for  the  elliptical  portion 
from  equation  (Al).  The  counterweight  lobes  are  then  treated  as  semi-circular  segments, 
and  their  contribution  to  mass  polar  moment  of  inertia  is: 


where  a  is  the  angle  subtended  by  the  counterweight  lobe;  120°  for  crankweb  (a)  and  70° 
for  crankweb  (b). 

Rotating  inertia  for  the  dynamometer  is  taken  from  Ref  [26].  The  coupling  shaft 
between  the  flywheel  and  the  dynamometer  is  calculated  as  a  circular  cylinder  using 
Equation  (Al). 

2.  Torsional  Rigidity 

For  modeling  of  the  torsional  rigidity,  the  components  of  the  crankshaft  must  be 
mathematically  converted  to  an  equivalent  shaft  of  a  constant  diameter.  This  is  a  simple 
geometric  problem  for  static  twisting  of  the  crankshaft,  but  becomes  very  complex  when 
considering  the  dynamic  crankshaft  twist  during  engine  operation. 

Wilson  [Ref  24]  presents  a  derivation  of  the  torsional  rigidity  for  various 
crankshaft  components.  The  rigidity  for  a  solid  cylindrical  shaft  of  diameter  D  is: 


K  = 


32  Le 


(A6) 


This  equation  can  be  used  to  determine  the  torsional  rigidity  of  any  solid  cylindrical 
component,  including  the  journal  bearings  and  crankpins.  For  the  crankweb,  an  effective 
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length  is  derived  based  on  the  bending  theory  of  beams.  The  crankweb  is  treated  as  a 
beam  subjected  to  a  bending  force  by  the  torsion  on  the  crankshaft: 


L . =  DA 


0.942 R 
T  W3 

,  iwhYrwh  . 


(A7) 


The  equivalent  shafting  between  concentrated  inertia  points  in  the  model  can  be 
determined  by  adding  up  the  torsional  rigidity  of  their  components.  However,  the 
equations  above  assume  an  unconstrained  crankshaft  deflection  due  to  an  applied  static 
torque.  The  true  equivalent  length  of  the  crankshaft  elements  will  be  modified  by  several 
other  factors:  local  deformation  where  the  journal  bearings  and  crankpins  join  the 
crankweb,  bearing  restraint  on  the  journal  bearing,  and  non-ideal  lever  arm  at  the 
crankweb  because  the  bearing  and  crankpin  are  not  attached  at  a  single  point.  [Ref  24] 

Two  empirical  relations  are  considered  in  this  study  to  account  for  increased 
rigidity  of  the  crankshaft  elements  due  to  dynamic  constrained  operation  while  mounted 
in  the  engine  block.  The  first  was  devised  by  Carter  [Ref  27]: 


4 ih=Di 


Ljh+0.$T 

K 


L  =D 

ecp 


0-75  Lq> 


D 


cp 


L  =D 

ecw 


L5R 
T  W 3 

1whYYwb 


(A8) 


and  the  second  by  Wilson  [Ref  24]: 


kjh  =  D4 


Ljh+0ADJh 

K 


L  =D 

cep 


L  +0.4  D 

cp  Cp 


D 


cp 


L  =D 

ecw 


R-0.2(D.fe  +  Dcp) 

T+WL 


(A9) 


Table  7  compares  the  results  obtained  by  each  of  the  methods  so  far  described  with  the 
values  from  the  model  supplied  by  Detroit  Diesel  [Ref  23].  The  general  assumptions  for 
the  two  empirical  methods  are  clear:  Carter’s  method  assumes  a  stiffer  crankpin  and  more 
flexible  crankwebs  while  Wilson’s  method  is  reversed  [Ref  24],  The  determination  of 
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torsional  rigidities  for  the  model  is  not  straightforward;  it  will  require  some  calibration 
from  measured  data. 


Table  7.  Calculated  Torsional  Rigidities 


Torsional  Rigidity  (106  lbf* in/rad)  by  method: 

Value 

Unconstrained 

Carter 

Wilson 

mUESlifl 

Kjb 

7.07 

41.5 

35.3 

KcD 

4.60 

38.4 

17.7 

Kcw 

31.5 

19.8 

49.4 

K, 

1.94 

3.11 

3.22 

3.47 

k2,  k3 

2.37 

6.61 

7.98 

9.55 

K4 

3.49 

10.8 

12.15 

13.50 

3.  Auxiliary  Loads 

The  auxiliary  loads,  with  the  exception  of  the  oil  pump,  are  driven  off  the  timing 
gear,  just  forward  of  the  flywheel.  For  the  model,  the  total  of  the  auxiliary  load  torque  is 
considered  to  be  placed  at  05.  The  contribution  of  the  individual  loads  can  be  determined 
by  deriving  the  power  required,  then  the  torque  is  related  to  the  power  P  by  the  equation: 


T  = 


2nN 


(A10) 


Two  cam  shafts  are  driven  off  the  timing  gear,  at  the  rear  of  the  crankshaft  and 
just  forward  of  the  flywheel.  The  load  due  to  these  camshafts  primarily  results  from  three 
components:  the  bearing  friction,  the  operation  of  the  fuel  injector  pistons,  and  the 
compression  of  springs  associated  with  the  injectors  and  the  valves.  Bearing  friction  is 
determined  in  the  next  section.  The  action  of  the  injector  pistons  is  considered  as  a 
polytropic  isothermal  process  of  compression  from  50  to  2800  psig.  Work  required  to 
compress  the  springs  is  calculated  from: 


=-K . 


spring 


spring 


(All) 
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where  5i  and  82  are  the  initial  and  final  spring  deflection,  respectively.  The  load  torque 
for  the  cam  shafts  is  calculated  from  the  work  done  by  the  shafts  over  one  rotation. 

The  fuel  pump  is  a  positive  displacement  gear  type  pump  that  provides  a  flowrate 
Q  of  1  gpm  at  65  psi  when  the  engine  is  at  2800  RPM.  Assuming  a  pressure  at  the  input 
of  about  15  psi,  the  torque  required  can  be  determined  from  the  power  P: 

P  =  Q&P  (A  12) 

Since  the  flowrate  Q  is  proportional  to  the  speed  N,  the  fuel  pump  torque  will  be  a 
constant  value  regardless  of  engine  speed. 

The  water  pump  is  a  centrifugal  pump  that  provides  37  gpm  at  2800  RPM.  Since 
a  pressure  drop  was  not  provided  from  the  service  manual,  its  effect  is  estimated  as 
comparable  to  the  other  auxiliary  components. 

A  roots  blower  provides  the  scavenging  pressure  that  clears  the  pistons  at  the 
bottom  of  the  stroke.  This  blower  is  rated  to  provide  338  cfm  at  2800  RPM. 

The  oil  pump  is  a  rotary  style  positive  displacement  pump  rated  to  deliver  15  gpm 
at  2800  RPM.  The  power  can  be  calculated  from 

P  =  mhA=(Qp){^A  (A  13) 

l  P  J 

where  hA  is  the  head  provided  by  the  pump  m  is  the  mass  flowrate  of  the  oil,  p  is  the  oil 
density,  and  Ap  is  the  pressure  increase.  The  increase  in  fluid  head  due  to  velocity 
increase  is  neglected. 

4.  Friction  Losses 

Journal  bearing  friction  is  accounted  for  by  assuming  that  the  two  surfaces  are 
completely  separated  by  the  lubricating  film;  that  hydrodynamic  lubrication  is  dominant. 
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Using  the  relations  provided  by  Heywood  [Ref  28],  an  equation  to  calculate  the  parasitic 
torque  due  to  a  bearing’s  friction  is  derived: 


2  Tfr  K2DhLiN  Wf 

Tf  =  fWfr  =>  /= — —  =  - — _£c__ ;  (j  —  — 1— 

f  J  WfDh  ho  LD, 


(A14) 


■b^h 


where  f  is  the  friction  factor,  Wf  is  the  bearing  load,  p  is  the  oil  viscosity,  h  is  the 
average  bearing  clearance,  N  is  the  rotation  speed  in  RPM,  and  Lb  and  Db  are  the  bearing 
length  and  diameter.  This  can  be  solved  for  torque  to  yield: 


Ah 


(A15) 


The  friction  torque  is  proportional  to  rotation  speed,  and  independent  of  bearing  load 
under  the  assumption  of  hydrodynamic  lubrication.  This  equation  can  be  applied  to  the 
crankshaft  main  bearings,  the  crankpins,  and  the  camshaft  bearings.  There  are  four  main 
bearings,  the  coefficients  Ci,  C2,  and  C3  in  the  equations  of  motion  are  determined  as 
one-third  of  the  total  friction  for  the  main  bearings.  Crankpin  friction  is  lumped  in  with 
the  piston  ring  friction,  and  the  friction  due  to  the  camshaft  bearings  is  included  with  the 
auxiliary  loads. 

Piston  ring  friction  is  not  easily  modeled  analytically,  but  is  instead  estimated  by 
empirical  methods.  From  Heywood  [Ref  28],  piston  friction  is  considered  as  the  sum  of 
two  components:  a  boundary  friction  generally  proportional  to  engine  loading,  and  a 
hydrodynamic  friction  proportional  to  piston  speed.  The  exact  relation  is  not  only 
difficult  to  predict  for  any  one  engine,  it  will  change  as  the  engine’s  condition  varies.  A 
more  detailed  study  of  the  ring  friction  is  beyond  the  scope  of  this  study.  Although  it  is 
expected  that  the  amount  of  piston  friction  varies  as  a  function  of  the  piston  speed,  a 
constant  parasitic  force  is  assumed,  which  corresponds  to  a  parasitic  torque  Tpar  which 
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varies  with  crank  angle,  as  in  Equation  (A  16).  The  derivation  of  the  force/torque 
relationship  is  detailed  in  the  next  Appendix.  Since  the  magnitude  of  Tpar  is  small 
relative  to  Tcyi,  this  is  considered  a  reasonable  approximation. 


T  —  F  R\ 

1  par  r  par1'' 


sin(0  +  (j) ) 


COS0 


(A  16) 


The  estimates  for  loads  and  frictions  formulated  above  are  meant  to  provide  a 
relative  relation  between  them.  For  this  study,  actual  losses  were  determined  from  the 
measured  pressure  data  (Table  8).  The  values  for  T|oad  and  Tpar  were  each  estimated  as  a 
fraction  of  the  total  losses. 


Table  8.  Empirical  Friction  Losses 


Engine 

Speed 

(RPM) 

Dyno 

Load 

(ft*lbf) 

Cyl  #1 
Average 
Gas 
Torque 

Cyl  #2 
Average 
Gas 
Torque 

Cyl  #3 
Average 
Gas 
Torque 

Total 

Average 

Gas 

Torque 

Total 

Losses 

Efficiency 

1000 

80 

24.9 

56.1 

132.5 

1000 

100 

32.4 

58.7 

62.5 

153.6 

53.6 

65.1  %  t 

mm 

73.1 

187.2 

52.2 

72.1  % 

mm 

79.0 

82.4 

215.1 

80.1 

62.8  % 

1500 

160 

63.2 

87.9 

91.9 

83.0 

65.8  % 

2000 

160 

72.7 

103.3 

118.6 

57.4  % 
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APPENDIX  B.  GEOMETRY  OF  RECIPROCATING  COMPONENTS 

The  nonlinear  motion  of  the  piston  and  connecting  rod  (Figure  21)  presents 
unique  geometrical  and  mathematical  problems  when  modeling  a  reciprocating  engine. 
The  following  derivations  are  common  in  the  literature,  but  are  presented  here  to 
document  the  detailed  analysis  required  to  formulate  the  torsional  model. 


Figure  21.  Piston  and  Connecting  Rod.  From  Ref  [19] 

1.  Indicated  Cylinder  Torque 

There  are  several  ways  to  derive  the  relation  between  the  gas  pressure  present  in 
the  cylinder  and  the  resulting  indicated  torque  applied  to  the  crankshaft.  Piston  ring 
friction  is  ignored  for  this  derivation;  it  is  corrected  for  by  a  parasitic  force  as  described 
in  the  previous  Appendix.  Figure  22  shows  the  relation  of  the  piston  to  the  crankshaft.  A 
static  analysis  assumes  two  forces  applied  at  the  piston  pin.  Fp  is  the  net  force  due  to  the 
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indicated  cylinder  pressure,  minus  the  reference  pressure  applied  to  the  underside  of  the 
piston  from  the  crankcase  (Fp  =  PnetAp).  Additionally,  a  reaction  force  Fr  is  applied  by  the 
cylinder  walls  on  the  piston  rings;  this  force  constrains  the  piston  to  linear  motion  in  the 


cylinder.  The  resultant  force  in  the  direction  of  the  connecting  rod  is: 


P  _  Ap  {Pcyi  Pre{  ) 

cr  “  COS0 


where  Ap  is  the  cross-sectional  area  of  the  piston,  Pcy]  is  the  indicated  cylinder  pressure, 
and  Pref  is  the  reference  pressure.  The  crank  angle  0  and  the  connecting  rod  angles  <|)  and 
y  are  related  by 

Lsin(<}>)  =  Rsin(0)  and  sin(y)  =  sin(0+<{>)  (B2) 

The  resultant  torque  applied  at  the  crankshaft  is  then  calculated  by  taking  the  cross- 
product  of  the  connecting  rod  resultant  force  vector  and  the  crankpin  position  vector,  and 
simplifying: 


T,=FxR=FjRsm 


r=Ar(p„,-P„,y— ^ 


COS0 


BDC  TDC 


Figure  22.  Geometry  of  Reciprocating  Components 
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The  same  result  is  also  achieved  by  equating  the  work  done  at  the  piston  with  the 
work  done  by  the  torque  on  the  crankshaft,  as  demonstrated  by  Taylor  [Ref  29]: 

P,„dV  =  l,,de  =>  T„,  --  P„,Ar  j  (B4) 

Taking  the  relation  for  s  from  Wilson  [Ref  24], 

s  =  Rcose  +  Lcos<p  =  Rcos0  +  [L2-R2  sin2  0\2  (B5) 

then 


_  .  _  R2sin6cos6 
-/(Sint/'-T-: - r - - — t 

=  6 

[-Rsin*-S!sin<>C0S',l 

(L2-R2  sin20j_ 

COS0 

and  using  trigonometric  relations,  equation  B4  will  simplify  to  equation  B3. 

2.  Reciprocating  Torque 

While  the  rotating  parts  of  the  crankshaft  maintain  a  nearly  constant  angular 
velocity,  the  reciprocating  components  are  alternately  accelerated  and  decelerated  in  a 
constrained  linear  motion.  At  the  crankshaft,  this  will  be  seen  as  a  load  torque  while  the 
piston  is  accelerated  from  TDC  to  its  maximum  speed,  and  will  supply  a  torque  as  the 
piston  is  decelerated  to  BDC.  Taylor  [Ref  29]  formulates  a  method  of  deriving  this 
reciprocating  torque. 

First,  it  is  necessary  to  determine  the  amount  of  the  reciprocating  mass.  Clearly, 
the  entire  mass  of  the  piston  contributes,  but  only  a  portion  of  the  connecting  rod  is 
reciprocating,  and  the  rest  must  be  considered  rotating  mass.  A  first  approximation 
idealizes  the  connecting  rod  as  two  lumped  masses  connected  by  a  massless  shaft,  with 
the  same  total  mass  and  center  of  gravity  as  the  real  connecting  rod  (Figure  23).  The 
center  of  gravity  for  the  real  connecting  rod  is  simply  found  by  balancing  the  rod;  for  this 
engine  h  =  3.5  in.  and  j  =  5.3  in.'  The  portion  labeled  Wi  is  then  added  into  the 
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reciprocating  mass,  while  the  remainder  is  considered  part  of  the  crankshaft  rotating 
mass.  Then,  the  instantaneous  torque  Trec  is  found  by  equating  the  change  in  the 
reciprocating  mass  kinetic  energy  with  the  work  done  by  the  crankshaft: 


(B7) 


Figure  23.  Idealized  Connecting  Rod.  From  Ref  [29] 

A  small  correction  must  then  be  made  to  account  for  the  difference  between  the 
polar  moment  of  inertia  for  the  idealized  connecting  rod  and  the  actual  polar  moment  of 
inertia: 


^ rec.corr  ^  cr 


w  i  y  RcosO 


g  )  Lcoscj) 


where  Jcr  is  the  polar  moment  of  inertia  for  the  actual  connecting  rod.  While 


Taylor  derives  series  relations  to  state  all  values  as  functions  of  0,  for  this  study  the  angle 
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§  and  the  distance  s  are  calculated  directly,  then  differentiated  numerically  within  the 
program. 

3.  Reciprocating  Inertia 

During  engine  operation,  the  reciprocating  components  contribute  to  the  rotating 
inertia  of  the  crankshaft.  However,  the  magnitude  of  the  reciprocating  inertia  varies  as  a 
function  of  the  crank  angle  0.  The  rotating  motion  of  the  crankshaft  drives  a  linear 
motion  of  the  piston  in  the  cylinder.  When  the  piston  is  at  TDC,  an  incremental  rotation 
of  the  crankshaft  results  in  zero  linear  motion  of  the  piston,  while  at  90°  the  same 
incremental  rotation  of  the  crankshaft  results  in  maximum  linear  motion  of  the  piston. 
Therefore,  the  influence  of  the  piston  mass  on  the  inertia  seen  at  the  crankshaft  will  vary 
during  crankshaft  rotation.  Normally,  crankshaft  inertial  models  include  an  average 
value  of  one-half  the  maximum  reciprocating  inertia  to  account  for  the  reciprocating 
components.  However,  in  this  application,  we  are  interested  in  crank-angle  dependent 
values  of  torque,  so  we  must  account  for  this  crank-angle  dependent  variation  of 
reciprocating  inertia.  The  reciprocating  inertia  can  be  calculated  as  a  function  of  the 
reciprocating  mass,  the  eccentricity,  and  the  crank  angle: 

W  R 2 

Jm.  =  (1  -  cos(20))  (B9) 

2  g 

The  value  Wrec  is  the  weight  of  the  reciprocating  components,  as  determined  previously. 

Reciprocating  inertia  is  a  separate  effect  from  the  reciprocating  torque  already 
described.  The  changing  value  of  reciprocating  inertia  accounts  for  the  extra  mass  that, 
along  with  the  rotating  mass,  must  be  accelerated  when  the  crankshaft  is  accelerated.  The 
reciprocating  torque  accounts  for  the  energy  required  to  accelerate  this  reciprocating 
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mass  from  zero  to  the  current  rotational  speed  of  the  crankshaft,  before  it  is  added  to  the 
rotating  mass. 
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APPENDIX  C.  NATURAL  FREQUENCY  AND  MODAL  ANALYSIS 

A  modal  analysis  can  be  conducted  for  the  crankshaft  using  the  torsional  model. 
Fourier  analysis  of  the  torsional  vibrations  measured  at  the  optical  encoder  will  show 
peaks  corresponding  to  the  measured  natural  frequencies.  Comparison  of  the  predicted 
natural  frequencies  from  the  model  to  these  measured  natural  frequencies  is  a  powerful 
calibration  tool  for  fine-tuning  the  model. 

The  matrix  equation  describing  the  torsional  mode  of  the  crankshaft  is  repeated 

here: 


[J]6+[C]0  +  [K]0  =  [T]  (7) 

Neglecting  the  damping  effects,  the  natural  frequencies  are  calculated  by: 

[*]- [j\o2  =  [0]=>  co1  =  eig{jy[K^  (Cl) 

For  the  given  model  parameters  (Table  1),  the  results  are  tabulated  in  Figure  24.  There 
are  six  modes  of  natural  vibration,  with  the  first  mode  being  the  trivial  rigid-body 
oscillation,  where  there  is  no  crankshaft  twist. 

A  natural  vibration  component  due  to  the  rigidity  and  inertia  of  the  flexible 
coupling  and  the  optical  encoder  disk  is  expected.  The  given  values  for  the  unit  are: 

Jrotor=  1-45  X  10'6  kg*m2 

J coupling  =  3  x  10'6  kg*m2 

Kcoupling  =  200  N*m/rad  [Ref  2 1  ] 

and  given 


the  natural  frequency  for  torsional  vibration  for  the  coupling/encoder  unit  is: 
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co  =  1067  Hz 


Not  surprisingly,  a  large  response  is  seen  in  the  measured  data  at  precisely  this  frequency. 


O.E.  #1  #2  #3  flywheel  dyno 


Figure  24.  Torsional  Vibration  Modes 


For  highest  resolution,  a  Fast  Fourier  Transform  (FIT)  is  conducted  on  the 
measured  angular  velocity  at  the  optical  encoder  over  the  entire  1 1  cycles  of  collected 
data.  Figure  25  shows  the  frequency  spectrum  for  measured  data  from  the  1000  RPM, 
100  ft*lbf  run.  Spectrums  obtained  for  other  engine  speeds  and  loads  are  similar.  Figure 
26  is  an  expanded  view  of  the  low-frequency  portion  of  the  spectrum.  A  “comb”  of 
amplitude  spikes  are  seen,  corresponding  to  harmonics  of  the  engine  rotation  frequency, 
16.7  Hz.  As  expected,  a  large  frequency  response  is  seen  at  about  1060  Hz, 
corresponding  to  the  natural  frequency  of  the  encoder/coupling  unit.  Additional 
responses  are  seen  at  about  430  Hz,  1 130  Hz,  and  1910  Hz,  and  these  agree  with  three  of 
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the  expected  modes  from  the  torsional  model.  A  peak  seen  at  about  2400  Hz  is  from  an 
unknown  source;  it  is  not  present  in  frequency  spectrums  taken  at  other  engine  speeds. 


Natural  frequencies  that  are  predicted  by  the  model  but  not  seen  in  the  spectrum 
are  probably  due  to  low  amplitudes  at  the  crankshaft  nose.  For  instance,  the  expected 
344  Hz  frequency  is  mostly  oscillation  of  the  dynamometer  rotor  with  respect  to  the 
flywheel;  the  crankshaft  oscillation  amplitude  would  be  much  smaller.  As  expected,  each 
of  the  predicted  modes  has  a  node  close  to  the  flywheel  because  it  contains  the  bulk  of 
the  system  inertia.  Measurements  taken  at  the  flywheel  would  be  expected  to  show 
almost  no  high  frequency  vibration,  and  this  is  seen  in  the  measured  data. 


Frequency 

Figure  25.  Frequency  Spectrum  for  Measured  Angular  Velocity  at  0i  (0-3000  Hz) 
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amplitude 


Crankshaft  nose  (Theta  one)  (rad)  Flywheel  (Theta  five)  (rad) 


APPENDIX  D.  ADDITIONAL  DATA  PLOTS 


x  1 0  3  1 000  RPM,  80  ft*!bf  Phase  Deviation 


4r 
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Time  (sec) 


xIO'3 
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Time  (sec) 

Figure  27.  Phase  Deviation  (1000  RPM,  80  Ft*Ibf) 
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Cyl#3  Torque  (ft*lbf)  Cyl  #2  Torque  (ft*lbf)  Cyl  #1  Torque  (ft*lbf) 


1500  RPM,  135  ft*lbf  Cylinder  Gas  Torques 
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Crank  Angle  (deg) 


Figure  39.  Individual  Cylinder  Gas  Torques  (1500  RPM,  135  Ft*Ibf) 
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Cyi#3  Torque  (ft*lbf)  cy’  *2  Torclue  (ft*lbf)  Cyl  #1  Torque  (ft 


1 500  RPM,  1 60  ft*lbf  Cylinder  Gas  Torques 


IiTiTil 


0  50  100  150  200  250  300  350 

Crank  Angle  (deg) 

Figure  40.  Individual  Cylinder  Gas  Torques  (1500  RPM,  160  Ft*lbf) 
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Total  Cylinder  Gas  Torque  (ft*lbf) 


1000  RPM,  80  ft*lbf  Total  Cylinder  Gas  Torque 


Crank  Angle  (deg) 


Figure  42.  Total  Gas  Torque  (1000  RPM,  80  Ft*lbf) 
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Total 


Total  Cylinder  Gas  Torque  (ft*lbf) 


1500  RPM,  160  ft*lbf  Total  Cylinder  Gas  Torque 


Crank  Angle  (deg) 


Figure  45.  Total  Gas  Torque  (1500  RPM,  160  Ft*lbf) 
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Total  Cylinder  Gas  Torque  (ft*lbf) 


Crank  Angle  (deg) 


Figure  46.  Total  Gas  Torque  (2000  RPM,  160  Ft*lbf) 
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APPENDIX  E.  MATLAB  CODES 


Program  used  to  calculate  cylinder  gas  torques  from  measured  pressures. 

%  TORQUE36 

%  This  code  computes  the  individual  torque  contribution  of  each 
%  individual  cylinder  referenced  to  TDC  of  Nr.l  Cylinder. 

%  Gas  torque  is  calculated  based  on  measured  pressures 

load  walklOO.md  %  ECA  cylinder  #1  pressure  data 

load  wblklOO.md  %  ECA  cylinder  #2  pressure  data 

load  wclklOO.md  %  ECA  cylinder  #3  pressure  data 

pa  =  reshape  (walkl00,5,720);  Plcyl  =  [pa(l ,:)  pa(l)]; 

pb  =  reshape  (wblkl00,5,720);  P2cyl  =  [pb(l ,:)  pb(l)]; 

pc  =  reshape  (wclkl00,5,720);  P3cyl  =  [pc(l ,:)  pc(l)]; 

W  =  7.556;  %  Reciprocating  weight  (lbf) 

R  =  2.25 ;  %  Crankshaft  Eccentricity  (in) 

B  =  3.875;  %  Cylinder  Bore  (in) 

L  =  8.8;  %  Connecting  Rod  length  (in) 

g  =  386;  %  Gravitational  acceleration  (lbf*in/secA2) 

theta  =  linspace(0,2*pi,721);  %  Crank  angle  vector 

N  =  1022;  %  RPM 

Load  =  100;  %  Ft*lbf 

omega  =  2*pi*N/60;  %  rad/sec 

dt  =  (60/N)/720; 

si  =  R*cos(theta)  +  sqrt(LA2  -  (RA2)*sin(theta).A2); 

s2  =  R*cos(theta-4*pi/3)  +  sqrt(LA2  -  (RA2)*sin(theta-4*pi/3)  A2); 

s3  =  R*cos(theta-2*pi/3)  +  sqrt(LA2  -  (RA2)*sin(theta-2*pi/3)  A2); 

Spl  =  deriv(sLdt);  Spdl  =  deriv(Spl,dt);%  Piston  speed  (in/sec)  and 
Sp2  =  deriv(s2,dt);  Spd2  =  deriv(Sp2,dt);%  Piston  acceleration  (in/secA2) 

Sp3  =  deriv(s3,dt);  Spd3  =  deriv(Sp3,dt); 

Tlrec  =  -(W/g)*Spdl  *(Spl/omega)/12; 

T2rec  =  -(W/g)*Spd2’.*(Sp2/omega)/12; 

T3rec  =  -(W/g)*Spd3.*(Sp3/omega)/12; 

Fpar  =  201 ;  %  lbf  Parasitic  force 

Tparl  =  Fpar*abs(Spl/omega)/12;  %  ft*!bf  Parasitic  torque 

Tpar2  =  Fpar*abs(Sp2/omega)/12; 

Tpar3  =  Fpar*abs(Sp3/omega)/12; 
pref  =  14.706  +  0.00205*(N-634);  %  psia 

Tlcyl=-((Plcyl.*pref-pref)*(pi*BA2/4).*Spl/omega)/12;  %  FT*LBF 

T2cyl=-((P2cyl.*pref-pref)*(pi*BA2/4).*Sp2/omega)/12;  %  FT*LBF 

T3cyl=-((P3cyl.*pref-pref)*(pi*BA2/4).*Sp3/omega)/12;  %  FT*LBF 

Tbrg  =  (0.0003403)*N;  %  FT-LBF  Bearing  friction  per  cylinder 

Taux  =  38;  %  FT-LBF  Valve  train  and  auxiliaries 

Tpmp  =  0.792;  %  FT-LBF  Oil  pump  load 

Ttot  =  Tlcyl+T2cyl+T3cyl-Tparl-Tpar2-Tpar3-3*Tbrg-Taux-Tpmp+Tlrec+T2rec+T3rec;%  FT-LBF 
figure  (1) 

degrees  =  theta*  180/pi; 

plot  (degrees, Plcyl, Tex’, degrees, P2cyl,,ko’,degrees>P3cyl,,k■,) 
axis([0, 360, 0,90]) 

title  (Cylinder  Pressures  referenced  to  TDC  of  #1  Cylinder5) 
xlabel  (Crank  Angle  (degrees)5) 
ylabel  (Cylinder  Indicated  Pressure  (bars,  absolute)5) 
legend  (Cyl#l  ’.CyM^CyMB’) 
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grid 

orient  landscape 
figure(2) 

plot(theta, TlcyfbX5, theta, T2cyl,  bO5,  theta,' T3cyl,5b.5) 
grid,  hold  on 

plot(theta,T  1  rec,  VX’,  theta, T2rec,  ’rO’,theta,T3rec,  ’r. ) 
plot  (theta, Ttot,  5g5, theta, Load  *ones(size(theta)), ’c  *) 

title(lndividual  Cylinder  Torque  input  Referenced  to  TDC  of  Nr.  1  Cylinder5) 

ylabelfGas  Torque  (FT-LB)5) 

xlabel(Degrees  after  TDC  of  NR.  1  Cylinder5) 

legend(5Cyl  #l’,5Cyl  #2’, 5Cyl  #30 

orient  landscape 

%  COMPUTE  AVERAGE  TORQUE  CONTRIBUTION  OF  EACH  CYLINDER 
T 1  cy  I_avg=trapz(theta,T  1  cy  I)/(2*pi) 

T2cyl_avg=trapz(theta,T2cyl)/(2*pi) 

T3cyl_avg=trapz(theta,T3cyl)/(2*pi) 

Avg_Torque_input=trapz(theta,Ttot)/(2*pi) 

Torque_out_to_TorqueJn  =  Load/Avg_Torque_input 


Programs  used  to  calculate  angular  positions,  given  cylinder  gas  torques 


(time-marching  method). 


%  MEASPRD3 

%  Concurrently  plots  measured  and  predicted  responses 

%  using  time-marching  direct  integration  method 

% 

%  Section  One  plots  the  measured  response 


load  wl  lklOO.md 
t  =  wllklOO(:,l); 
tt  =  diff(t); 

tt  =  (reshape(tt, 720,1 1))5; 
tttt=  mean(tt,l); 

position  =  linspace(0,(2*pi),721); 
pos  =  position(l:720); 
omega  =  (l/sum(tttt))*2*pi; 
timem  -  [0  cumsum(tttt)]; 
angposm  =  position-timem.*omega; 
plot(timem,angposm,5rO’) 
hold  on 


%  Load  optical  encoder  data  file 

%  determine  dt’s 
%  Phase  lock  one  cycle 
%  Ensemble  average  the  phases 
%  The  known  positions  of  the  O.E.  windows 

%  Mean  rotational  velocity  (rad/sec) 

%  Time  vector  corresponding  to  position 
%  Angular  position  (radians) 

%  Plot  measured  angular  position  vs  time 


%  Section  Two  plots  the  predicted  response  based  on  pressure  data 
%  Uses  a  SIX  SECOND  ORDER  SUMULTANEOUS  EQUATION  ODE  SOLVER 
%  Calls  "deqns.m"  which  defines  the  system  of  2nd  order  ode’s 

% 

global  Tload  Tlcyl  T2cyl  T3cyl  Taux  Tpmp  j2rec  j3rec  j4rec  i 
load  walklOO.md  %  ECA  cylinder  #1  pressure  data 

load  wb lklOO.md  %  ECA  cylinder  #2  pressure  data 

load  wc lklOO.md  %  ECA  cylinder  #3  pressure  data 

pa  =  reshape  (walkl00,5,720);  pa  =  pa(l,:); 
pb  =  reshape  (wblkl00,5,720);  pb  =  pb(l ,:); 
pc  =  reshape  (wclkl00,5,720);  pc  =  pc(l,:); 

Plcyl  =  [pa(l,:)  pa(l)];%  Cylinder  pressures  (bars) 
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P2cyl  =  [pb(  1 ,:)  pb(l)]; 

P3cyl  =  [pc(l,:)  pc(  1 )] ; 
shp  =  [0  0  000-1. 8]*le-03; 
shv  =  [11111  1]*106.3; 
ic  =  [shp  shv]; 

Tload=1200; 

timep=Iinspace(0,sum(tttt),72I); 
angposp  =  zeros(l,721); 

B  =  3.875; 

R  =  2.25; 

W  =  7.556; 
g  =  386; 

L  =  8.8; 


%  define  IC’s 
%  lbf*inLoad  torque 
%  divide  one  rev  into  720  divisions 
%  initialize  predicted  position  vector 
%  in  Cylinder  bore 

%  in  Crankshaft  eccentricity 

%  lbf  Weight  of  reciprocating  components 

%  in/secA2  gravitational  acceleration 
%  in  Connecting  rod  length 


)  Piston  acceleration  (in/secA2) 


N  =  omega*60/(2*pi);  %  RPM 

pos  =  linspace  (0,(2*pi),721); 
pref  =  14.706  +  0.00205  *(N-634);  %  psia 
dt  =  (2*pi/omega)/720; 

si  =  R*cos(pos)  +  sqrt(LA2  -  (RA2)*sin(pos).A2); 
s2  =  R*cos(pos-4*pi/3)  +  sqrt(LA2  -  (RA2)*sin(pos-4*pi/3)  A2); 
s3  =  R*cos(pos-2*pi/3)  +  sqrt(LA2  -  (RA2)*sin(pos-2*pi/3).A2); 

Spl  =  derives l,dt);  Spdl  =  deriv(Spl,dt);  %  Piston  speed  (in/sec)  and 
Sp2  =  deriv(s2,dt);  Spd2  =  deriv(Sp2,dt); 

Sp3  =  deriv(s3,dt);  Spd3  =  deriv(Sp3,dt); 

Fpar=100;  %  lbf*in  Parasitic  force 

Tparl  =Fpar*abs(Spl /omega);  %  lbf*in  Parasitic  torque 

Tpar2  =  Fpar*abs(Sp2/omega); 

Tpar3  =  Fpar*abs(Sp3/omega); 

pos  =  linspace  (0,(2*pi),721); 

pref  =  14.706  +  0.00205 *(N-634);  %  psia 

Tpmp  =  9.5;  %lbf*in  Oil  pump  torque 

Tlcyl  =  -(Plcyl.*pref-pref)*(pi*BA2/4).*Spl/omega;  %  (in*lbf) 

T2cyl  =  -(P2cyl.*pref-pref)*(pi*BA2/4).*Sp2/omega;  %  (in*lbf) 

T3cyl  =  -(P3cyl.*pref-pref)*(pi*BA2/4).*Sp3/omega;  %  (in*lbf) 

Taux  =  168;  %lbf*in  Valvetrain  and  auxilliary  torque 
j2rec  =  (W*RA2/(2*g))*(l-cos(2*pos));  %lbf*in*secA2 

j3rec  =  (W*RA2/(2*g))*(l-cos(2*(pos-4*pi/3)));  %lbf*in*secA2 

j4rec  =  (W*RA2/(2*g))*(l-cos(2*(pos-2*pi/3)));  %lbf*in*secA2 

Tlrec  =  -(W/g)*Spd  1 .  *(Sp  1/omega); 

T2rec  =  -(W/g)*Spd2.*(Sp2/omega); 

T3rec  =  -(W/g)*Spd3.*(Sp3/omega); 

Tlcyl  =  Tlcyl  -  Tparl  +  Tlrec; 

T2cyl  =  T2cyl  -  Tpar2  +  T2rec; 

T3cyl  =  T3cyl  -  Tpar3  +  T3rec; 
step  =  1 ; 

theta  =  zeros(721,6); 
theta(l,:)  =  shp; 
for  i  =  l:step:720; 

[T,x]  =  ode45(’deqns’,[timep(i)  timep(i+step)],ic); 

ic  =  x(length(T),:);  %  reinitialize  IC’s  from  previous  iteration 

theta(i+step,:)  =  x(length(T),l:6); 
angposp(i+step)  =  x(length(T),l)  -  omega*timep(i+step); 
end 

plot(timep, angposp, ’bX’)  %  Plot  predicted  angular  position  vs  time 

title(Time  Marching  ODE45  Method  Comparison  to  Measured  Data*) 
xlabel(Time  (sec)’) 
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ylabeI(Theta  one  Phase  Deviation  (radians)’) 

legend  (Meas  ’,Pred  ) 

grid 

orient  landscape 
shaft_pos  =  (ic(l:6)-2*pi)’ 
cycle_stat  =  (ic  -  [(shp+2*pi)  shv])’ 

%  Plot  predicted  angular  velocity 
figure  (2) 

omegap  =  [shv(l)  (diff(angposp)./timep(2)+omega)]; 

omegam  =  (l./(720*tttt))*2*pi; 

omegam  =  [omegam  omegam(720)]; 

plot(timep, omegap,  ,b’,timep, omegam,  ’r) 

hold  on 

plot(timep,omega*ones(size(timep))) 
xlabel(’time(sec)) 
ylabeI(’Angular  Velocity  (rad/sec)) 
grid,  orient  tall 
figure  (3) 

plot(timem,angposm,  Tc. ) 
hold  on 

p!ot(timep,angposp,’k) 

titie(Time  Marching  ODE45  Method  Comparison  to  Measured  Data7) 
xlabel(Time  (sec)) 

ylabel(Theta  one  Phase  Deviation  (radians)) 

legend  (Meas  ’,Pred  ) 

grid 


%  DEQNS  function  to  determine  six  second  order  differential  equations 
%  To  be  used  in  ode45  fctn  in  measpred  program 


function  xdot=deqns(t,x) 

% - 

global  Tload  Tlcyl  T2cyl  T3cyl  Taux  Tpmp  j2rec  j3rec  j4rec  i 
% - Constants  to  be  used  in  differential  equations - 


j  1=02443; 
j2=2482; 
j3=.  1462; 
j4=.2482; 
j5=7.222; 
j6=.2870; 
kl=3.1  le6; 
k2=7.00e6; 
k3=7.00e6; 
k4=10.82e6; 
k5=1.304e6; 
cl2=.01 ; 
c23=.01; 
c34=01; 
c45=.01 ; 
c56=.01; 
c2=0.013; 
c3=0.013; 
c4=0.013; 

% - 

%  1 2  first  order  equations  which  define  the 


%Ib*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 
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%  original  6  second  order  equations  of  motion 

% - 

xdo(l)=x(7); 

xdo(2)=x(8); 

xdo(3)=x(9); 

xdo(4)=x(10); 

xdo(5)=x(ll); 

xdo(6)=x(12); 

xdo(7)=((c  1 2*x(8))-(c  1 2*x(7))+(k  1  *x(2))-(kl  *x(  1  ))-Tpmp)/j  1 ; 
xdo(8)=(Tlcyl(i)+(c23*x(9))-((c23+cl2+c2)*x(8))+(cl2*x(7))+(k2*x(3))... 

-((kl+k2)*x(2))+(kl*x(l)))/(j2+j2rec(i)); 

xdo(9)=(T2cyI(i)+(c34*x(10))-((c34+c23+c3)*x(9))+(c23*x(8))+(k3*x(4))... 

-((k2+k3)*x(3))+(k2*x(2)))/(j3+j3rec(i)); 
xdo(  1 0)=(T3cyl(i)+(c45  *x(  1 1  ))-((c45+c34+c4)*x(  1 0))+(c34*x(9))+(k4  *x(5)).  .. 

-((k3+k4)*x(4))+(k3*x(3)))/(j4+j4rec(i)); 
xdo(  1 1  )=((c56*x(  1 2))+((c56+c45)*x(  1 1  ))+(c45*x(  1 0))+(k5  *x(6))..; 

-((k4+k5)*x(5))+(k4*x(4))-Taux)/j5; 
xdo(  1 2)=(-Tload-(c56*x(  1 2))+(c56*x(  1 1  ))-(k5  *x(6))+(k5*x(5)))/j6; 
xdot=xdo’;  %  vector  defining  the  equations  of  motion 

Programs  used  to  calculate  angular  positions,  given  measured  cylinder  gas 
torques  (finite  element  method). 

%  MEASPS 

%  Determines  and  plots  comparisons  of  the  following: 

%  (1)  Measured  response  from  flywheel  and  optical  encoder  data 

%  (2)  Predicted  response  from  inertial  model  based  on  measured 

%  pressure  data  from  the  three  cylinders 

%  Evaluates  predicted  response  using  a  finite  element  formulation 

% - 

%  Section  One  plots  the  measured  response 

% - 

% - Measured  response  from  flywheel - 

load  w5 lklOO.md  %  Load  time  data  from  MDA 

t=  [0;w51kl00(l:7938,l)];  %  Extract  time  data  only 

tt  =  diff(t);  %  COMPUTES  THE  DELJTS 

tt  =  (reshape(tt,  1 26,63))’;  %  Phase  lock  one  cycle 

tttt  =  mean(tt);  %  Ensemble  average  the  phases 

dtrat  =  tt(l)/mean(tt(2:63,l)); 

teeth  =  linspace(0,2*pi,127);  tooth  =  2*pi/127; 

th51  =  tooth  *dtrat; 

pos5  =  [0  teeth(  1 : 1 25  )+th5 1  ] ; 

pos5(127)  =  2*pi; 

time5  =  [0,cumsum(tttt)-(l  -dtrat)  *tttt(  1 )] ; 
time5(127)  =  sum(tttt); 
omega  =  (l/max(time5))*2*pi; 

angpos5  =  pos5-time5.*omega;%Angular  position  (radians) 
figure  (1) 
subplot  (2,1,1) 

plot(time5,angpos5,’k.’)  %  Plot  measured  angular  position  vs  time 

hold  on 

% - Measured  response  from  optical  encoder - 

load  wl  lklOO.md  %  Load  time  data  from  MDA 

t  =  [0;wl  lkl00(l:7920,l)];  %  Extract  time  data  only 
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tt  =  diff(t); 

tt  =  (reshape(tt,720,ll))’; 
tttt  =  mean(tt); 

posl  =  Iinspace(0,(2*pi),721); 

timel  =  [0  cumsum(tttt)] ; 

timel  =  timel  *(max(time5)/max(timel)); 

angposl  =  posl -timel.  *omega; 

figure  (1) 

subplot  (2,1,2) 

plot(time  1  ,angpos  1 , k.  *) 

hold  on 


%  COMPUTES  THE  DEL.T’S 
%  Phase  lock  one  cycle 
%  Ensemble  average  the  phases 
%  Crank  angle  position  for  O.E.  windows 
%  Time  vector  corresponding  to  position 
%  Adjust  times  to  agree 
%  Angular  position  (radians) 

%  Plot  measured  angular  position  vs  time 


%  Section  Two  plots  the  predicted  response  based  on  pressure  data 

% - 

global  j2rec  j3rec  j4rec 

% . Cylinder  pressures  (bars) - 

load  walklOO.md  %  EC  A  cylinder  #1  pressure  data 

load  wblklOO.md  %  ECA  cylinder  #2  pressure  data 

load  wclklOO.md  %  ECA  cylinder  #3  pressure  data 

pa  =  reshape  (walklOO, 5,720);  Plcyl  =  [pa(l,:)  pa(l)]; 
pb  =  reshape  (wblkl00,5,720);  P2cyl  =  [pb(l ,:>  pb(l)]; 
pc  =  reshape  (wc lk 100, 5, 720);  P3cyl  =  [pc(l ,:)  pc(l)]; 

%  Variable  descriptions 
%  k  -  element  matrix 
%  f  =  element  vector 
%  kk  =  compressed  system  matrix 
%  ff  =  system  vector 

%  bcdof  =  a  vector  containing  dofs  associated  with  boundary  conditions 
%  bcval  =  a  vector  containing  boundary  condition  values  associated  with 
%  the  dofs  in  ’bcdof 


input  data  for  control  parameters 


nel  =  720; 

%  number  of  elements 

nnel  =  2; 

%  number  of  nodes  per  element 

ndof  =  6; 

%  number  of  dofs  per  node 

nnode  =  72 1 ; 

%  total  number  of  nodes  in  system 

sdof  =  nnode*ndof; 

%  total  system  dofs 

Tload  =  1200; 

%  (lbf*in) 

R  =  2.25; 

%  (in)  Crankshaft  eccentricity 

W  =  7.556; 

%  (lbf)  Weight  of  reciprocating  components 

g  =  386; 

%  (in/secA2)  Gravitational  acceleration 

B  =  3.875; 

%  (in)  Cylinder  Bore 

L  =  8.8; 

%  (in)  Connecting  Rod  length 

N  =  omega*60/(2*pi); 
dt  =  (2*pi/omega)/720; 

%  RPM 

si  =  R*cos(posl)  +  sqrt(LA2  -  (RA2)*sin(posl).A2); 

s2  =  R*cos(posl-4*pi/3)  +  sqrt(LA2  -  (RA2)*sin(posl-4*pi/3).A2); 

s3  =  R*cos(posl-2*pi/3)  +  sqrt(LA2  -  (RA2)*sin(posl-2*pi/3).A2); 

Spl  =  deriv(sl,dt);  Spdl  =  defiv(Spl,dt);  %  Piston  speed  (in/sec)  and 

Sp2  =  deriv(s2,dt);  Spd2  =  deriv(Sp2,dt);  %  Piston  acceleration  (in/secA2) 

Sp3  =  deriv(s3,dt);  Spd3  =  deriv(Sp3,dt); 

Fpar  =  100; 

%  lbf*in  Parasitic  force 

Tparl  =  Fpar*abs(Sp  1/omega); 
Tpar2  =  Fpar*abs(Sp2/omega); 

%  lbf*in  Parasitic  torque 
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Tpar3  =  Fpar*abs(Sp3/omega); 
pos  =  linspace  (0,(2*pi),721); 
pref  =  14.706  +  0.00205*(N-634);  %  psia 
Tpmp  =  9.5;  %lbf*in  Oil  pump  torque 
Tlcyl  =  -(Plcyl.*pref-pref)*(pi*BA2/4).*Spl/omega; 

T2cyl  =  -(P2cyl.*pref-pref)*(pi*BA2/4).*Sp2/omega; 

T3cyl  =  -(P3cyl.*pref-pref)*(pi*BA2/4).*Sp3/omega; 

Taux  =  168;  %lbf*in  Valvetrain  and  auxiliary  torque 
j2rec  =  (W*RA2/(2*g))*(l-cos(2*pos)); 
j3rec  =  (W*RA2/(2*g))*(l-cos(2*(pos-4*pi/3))); 
j4rec  =  (W*RA2/(2*g))*(l-cos(2*(pos-2*pi/3))); 

Tlrec  =  -(W/g)*Spdl.*(Spl/omega); 

T2rec  =  -(W/g)*Spd2.*(Sp2/omega); 

T3rec  =  -(W/g)*Spd3.*(Sp3/omega); 

Tlcyl  =  Tlcyl  -  Tparl  +  Tlrec; 

T2cyl  =  T2cyl  -,Tpar2  +  T2rec; 

T3cyl  =  T3cyl  -  Tpar3  +  T3rec; 
pack 

% - 

%  input  data  for  nodal  coordinate  values 

% - 

tcoord  =  linspace(0,max(time5),721); 

% - 

%  input  data  for  nodal  connectivity  for  each  element 

% - 

nodes  =  [(l:nel)’,(2:nnode)]; 

% - 

%  input  data  for  boundary  conditions 

% - 

%  Dirichlet  Boundary  Conditions 
shp  =[0  0  0  0  0  -200]*le-05; 

bcdof  =  [l,2,3,4,5,6,sdof-5,sdof-4,sdof-3,sdof-2,sdof-l,sdof]; 
bcval  =  [shp  (shp+2*pi)]; 

% - 

%  initialization  of  matrices  and  vectors 

% - 

ff  =  zeros(sdof,l);  %  initialization  of  system  force  vector 

kk  =  zeros(sdof,  1 8);  %  initialization  of  compressed  system  matrix 

index  =  zeros(sdof.l);  %  initialization  of  kk  index  vector 

% - 

%  computation  of  element  matrices  and  vectors  and  their  assembly 

% - 

for  i=l:nel;  %  loop  for  the  total  number  of  elements 

nl=nodes(i,l);  nr=nodes(i,2);  %  extract  nodes  for  (iel)-th  element 
tl=tcoord(nl);  tr=tcoord(nr);  %  extract  nodal  coord  values  for  the  element 
k  =  Kelm(tl,tr,i);  %  compute  element  matrix 

f  =  ((tr-tl)/2) * [-T pmp;  Tlcyl(i);  T2cyl(i);  T3cyl(i);  -Taux;  -Tload;... 

-Tpmp;  Tlcyl(i+1);  T2cyl(i+1);  T3cyl(i+1);  -Taux;  -Tload]; 

%  compute  element  vector 

for  ii  =  1 : 12;  %  assemble  element  matrices  and  vectors 

ff((i-l)*6+ii)  =  ff((i-l)*6+ii)  +  f(ii); 
if  ii<=6; 

kk((i-l)*6+ii,:)  =  kk((i-l)*6+ii,:)  +  [0,0, 0,0,0, 0,k(ii,:)]; 
index(i*6+ii)  =  (i-l)*6; 
else 


%  (in*lbf) 
%  (in*lbf) 
%  (in*lbf) 


%lbf*in*secA2 

%lbf*in*secA2 

%lbf*in*secA2 
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kk((i-l)*6+ii,:)  =  kk((i-l)*6+ii,:)  +  [k(ii,:),0, 0,0,0, 0,0]; 
end 
end 
end 

index(l:6)  =  -6*[l;l;l;l;l;l]; 

% - 

%  apply  boundary  conditions 

% - 

%  Dirichlet  Boundary  Conditions 
for  i  =  l:length(bcdof); 

kk(bcdof(i),:)  =  zeros(l,18); 
kk(bcdof(i),bcdof(i)-index(bcdof(i)))  =  1; 
ff(bcdof(i))  =  bcval(i); 
end 


%  Transform  kk  and  ff  to  upper  diagonal 

% - 

for  i  =  l:(sdof-l); 
f=fix((i-l)/6)*6+12; 
if  i>(sdof-6);  f=sdof;  end 

for  ii  -  (i+l):f;  %  other  rows  with  data  in  column  i 

v  =  kk(ii,i>index(ii))./kk(i,i-index(i));  %  multiplier 
kk(ii,i-index(ii))  =  0; 

for  j  =  (i+l):f;  %  data  elements  in  row  ii 

kk(ii  j-index(ii))  =  kk(ii  j-index(ii))-v*kk(i,j-index(i)); 
end 

ff(ii)  =  ff(ii)  -  v*ff(i); 
end 
end 


%  Solve  matrix  eqn  for  theta 

% - 

theta  =  zeros(l,sdof); 

theta(sdof)  =  ff(sdof)/kk(sdof,sdof-index(sdof)); 
for  i  =  (sdof-l):(-l):l; 
f  =  fix((i-l)/6)*6+12; 
if  i>(sdof-6);  f=sdof;  end 
for  j  =  (i+l):f; 

ff(i)  =  ff(i)  -  kk(i,j-index(i))*theta(j); 
end 

theta(i)  =  ff(i)/kk(i,i-index(i)); 
end 

theta  =  reshape(theta,6,nnode); 
angpospl  =  theta(l,:)  -  omega*tcoord; 
angpospS  =  theta(5>:)  -  omega*tcoord; 

% . . . 

%  Plot  predicted  angular  position  vs  time 

% . . . 

figure  (1) 
subplot(2,l,l) 
plotCtcoord^ngpospS^) 
xlabeI(Time  (sec)) 

title  (’1000  RPM,  100  ft*lbf  Phase  Deviation) 

ylabel(Tlywheel  (Theta  five)  (rad)) 

grid 
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subp!ot(2,l,2) 
plot(tcoord,angpospl  ,1c7) 
xlabel(Time  (sec)7) 

ylabel( ’Crankshaft  nose  (Theta  one)  (rad)7) 
grid 

legend  (Meas  ’,Pred  *) 

% - 

%  Compare  Other  degrees  of  freedom  to  theta  5 

% - - - 

twist  =  theta  -  ([  1 ;  1 ;  1 ;  1 ;  1 ;  1  ] *theta(5,:)); 
figure  (3) 

plot  (tcoord,twist(l,:),’y’),  hold  on 
plot  (tcoord,twist(2,:),’rT) 
plot  (tcoord,twist(3,:),’g’) 
plot  (tcoord,twist(4,:),’b’) 
plot  (tcoord,twist(6,:),’cT) 

title( ’Angular  Deviation  from  Theta  Five  (Flywheel)7) 
xlabel(Time  (sec)7) 

ylabel(’Angular  Deviation  from  Theta  Five  (radians)7) 
legend  (’O.E.VCyl  #1 ’,’Cyl  #2’,’Cyl  #3’, Dyno’ Measured7) 
grid 

orient  landscape 

% - 

%  Plot  Measured  vs.  Predicted  Crankshaft  Twist 

% - 

thetarl  =  interpl  (time  l,posl,tcoord, ’spline7); 
thetar5  =  interpl  (time5, pos5,tcoord,  ’spline7); 
figure  (6) 

plot  (tcoord,thetarl-thetar5,’k.’) 

grid,  hold  on 

plot  (tcoord^wis^,:),^7) 

title  (’1000  RPM  100  ft*lbf  Crankshaft  Twist7) 

xlabel  (Time  (sec)7) 

ylabel  (Twist  (radians)7) 

legend(Meas  ’,Pred  0 

function  [k]  =  kelm(tl,tr,i) 

% 

%  KELM  calculates  the  element  matrix  k, 

%  given  tl  and  tr  as  inputs 
%  Global  variables 
global  j2rec  j3rec  j4rec 

% - 

%  Constants 

% - 

h  =  tr-tl; 
jl  =0.02443; 
j2  =  0.2482; 
j3  =0.1462; 
j4  =  0.2482; 
j5  =  7.222; 
j6  =  0.2870; 
kl  =  3.1  le6; 
k2  =  7.00e6; 
k3  =  7.00e6; 


%lbf*in*secA2/rad 
%lbf  *in  *sec  A2/rad 
%lbf*in*secA2/rad 
%lbf*in*secA2/rad 

%lbf*in*secA2/rad 
%lbf*in*secA2/rad 
%lbf*in/rad 
%lbf*in/rad 
%lbf*  in/rad 
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k4  =  10.82e6; 
k5=  1.304e6; 
cl2  =  0.01; 
c23  =  0.01; 
c34  =  0.01; 
c45  =  0.01; 
c56  =  0.01; 
c2  =  0.013; 
c3  =  0.013; 
c4  =  0.013; 

% - 


%lbf*in/rad 

%lbf*in/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 

%lbf*in*sec/rad 


k  =  zeros  (12,12); 


k(  1 , 1  )=  l/6*(2*kl  *trA3-2*k  1  *tl  A3-3*c  1 2*trA2-3  *c  1 2*tlA2+6*kl  *tr*tlA2-6*j  1  *tr+6*. . . 
j  1  *tl+6*c  1 2*tr*tl-6*kl*trA2*tl)/hA2; 

k(l,2)=l/6*(-2*kl*trA3+2*kl*tlA3+3*cl2*trA2+3*cl2*tlA2-6*kl*tr*tlA2-6*cl2*tr*tl+. 

6*kl*trA2*tl)/hA2; 

k(l,7)=-l/6*(-kl*trA3+kl*tlA3-3*cl2*trA2-3*cl2*tlA2-3*kl*tr*tlA2+3*kl*trA2*tl-6*... 
j  1  *tr+6*j  1  *tl+6*cl2*tr*tI)/hA2; 

k(l,8)=-l/6*(kl*trA3-kl*tIA3+3*cl2*trA2+3*cl2*tlA2+3*kl*tr*tlA2-3*kl*trA2*tl-6*... 

cl2*tr*tI)/hA2; 

k(2,l)=l/6*(-2*kl*trA3+2*kl*tlA3+3*cl2*trA2+3*cl2*tlA2-6*kl*tr*tlA2-6*cl2*tr*tl+. 

6*kl*trA2*tI)/hA2; 

k(2,2)=l/6*(6*k2*tr*tlA2-3*c2*trA2-3*c2*tlA2-6*j2rec(i)*tr+6*j2rec(i)*tl-6*j2*tr+... 

6*j2*tl+6*tl*tr*c23+6*tr*c2*tl-6*k2*tl*trA2+6*cl2*tr*tl-6*kl*trA2*tl+2*... 

k2*trA3-2*k2*tlA3-3*c23*trA2-3*c23*tlA2+2*kl*trA3-2*kl*tlA3-3*cl2*trA2-... 

3  *c  1 2*tl  A2+6  *k  1  *tr*tl  A2)/h  A2 ; 

k(2,3)=l/6*(-2*k2*trA3+2*k2*tlA3-6*k2*tr*tlA2+3*c23*trA2+3*c23*tlA2-6*tl*tr*c23+. 

6*k2*tI*trA2)/hA2; 

k(2,7)=-l/6*(kl*trA3-kl*tlA3+3*cl2*trA2+3*cl2*tlA2+3*kl*tr*tlA2-3*kl*trA2*tl-6*... 

cl2*tr*tl)/hA2; 

k(2,8)=-l/6*(-3*k2*tr*tIA2-3*c2*trA2-3*c2*tlA2-6*j2rec(i)*tr+6*j2rec(i)*tl-6*j2*... 

tr+6*j2*tl+6*tl*tr*c23+6*tr*c2*tl+3*k2*tl*trA2+6*cl2*tr*tl+3*kl*trA2*... 

tl-k2*trA3+k2*tlA3-3*c23*trA2-3*c23*tlA2-kl*trA3+kl*tlA3-3*cl2*trA2-3*... 

cl2*tlA2-3*kl*tr*tlA2)/hA2; 

k(2,9)=-l/6*(k2*trA3-k2*tlA3+3*k2*tr*tlA2+3*c23*trA2+3*c23*tlA2-3*k2*tl*trA2-6*... 

tl*tr*c23)/hA2; 

k(3,2)=l/6*(-2*k2*trA3+2*k2*tlA3-6*k2*tr*tlA2+3*c23*trA2+3*c23*tlA2-6*tl*tr*c23+. 

6*k2*tl*trA2)/hA2; 

k(3,3)=l/6*(6*k2*tr*tlA2+6*tl*tr*c23-6*k2*tl*trA2+6*k3*tr*tlA2-6*k3*trA2*tl+6*... 

tr*tl*c34+2*k2nrA3-2*k2*tlA3-3*c23*trA2-3*c23*tlA2+2*k3*trA3-2*k3*... 

tlA3-3*c34*trA2-3*c34*tIA2-3*c3*trA2-3*c3*tlA2+6*tl*c3*tr-6*j3*tr+6*... 

j3*tl-6*j3rec(i)*tr+6*j3rec(i)*tl)/hA2; 

k(3,4)=l/6*(-2*k3*trA3+2*k3*tlA3-6*k3*tr*tlA2+3*c34*trA2+3*c34*tlA2-6*tr*tl*c34+. 
6*k3  *tr  A2  *tl  )/hA2 ; 

k(3,8)=-l/6*(k2*trA3-k2*tlA3+3*k2*tr*tlA2+3*c23*trA2+3*c23*tlA2-3*k2*tl*trA2-6*... 

tl*tr*c23)/hA2; 

k(3,9)=-l/6*(-3*k2*tr*tlA2+6*tl*tr*c23+3*k2*tl*trA2-3*k3*tr*tlA2+3*k3*trA2*tl+6*... 

tr*tl*c34-k2*trA3+k2*tlA3-3*c23*trA2-3*c23*tlA2-k3*trA3+k3*tlA3-3*c34*... 

trA2-3*c34*tlA2-3*c3*trA2-3*c3*tlA2+6*tl*c3*tr-6*j3*tr+6*j3*tl-6*... 

j3rec(i)*tr+6*j3rec(i)*tl)/hA2; 

k(3,10)=-l/6*(k3*trA3-k3*tlA3+3*k3*tr*tlA2+3*c34*trA2+3*c34*tlA2-3*k3*trA2*tl-... 

6*tr*tl*c34)/hA2; 

k(4,3)=l/6*(-2*k3*trA3+2*k3*tlA3-6*k3*tr*tlA2+3*c34*trA2+3*c34*tlA2-6*tr*tl*c34+. 

6*k3*trA2*tl)/hA2; 
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k(4,4)=l/6*(6*tr*c4*tl+6*tr*c45*tl+6*k3*tr*tlA2-6*k3*trA2*tl+6*tr*tl*c34+2*k3*... 

trA3-2*k3*tlA3-3*c34*trA2-3*c34*tlA2+6*tr*k4*tlA2-6*k4*tl*trA2+2*k4*... 

trA3-2*k4*tlA3-3*c45*trA2-3*c45*tlA2-3*c4*trA2-3*c4*tlA2-6*j4*tr+6*... 

j4*tl-6*j4rec(i)*tr+6*j4rec(i)*tl)/hA2; 

k(4,5)=l/6*(-2*k4*trA3+2*k4*tlA3-6*tr*k4*tlA2+3*c45*trA2+3*c45*tlA2-6*tr*c45*... 

tl+6*k4*tl*trA2)/hA2; 

k(4,9)=-l/6*(k3*trA3-k3*tlA3+3*k3*tr*tlA2+3*c34*trA2+3*c34*tlA2-3*k3*trA2*tl-6*... 

tr*tl*c34)/hA2; 

k(4,10)=-l/6*(6*tr*c4*tl+6*tr*c45*tI-3*k3*tr*tlA2+3*k3*trA2*tl+6*tr*tl*c34-k3*... 

trA3+k3*tlA3-3*c34*trA2-3*c34*tlA2-3*tr*k4*tlA2+3*k4*tl*trA2-k4*trA3+... 

k4*tlA3-3*c45*trA2-3*c45*tlA2-3*c4*trA2-3*c4*tlA2-6*j4*tr+6*j4*tl-6*... 

j4rec(i)*tr+6*j4rec(i)*tl)/hA2; 

k(4,ll)=-l/6*(k4*trA3-k4*tlA3-3*k4*tl*trA2+3*c45*trA2+3*c45*tlA2+3*tri|tk4*tlA2-6*... 

tr*c45*tl)/hA2; 

k(5,4)=l/6*(-2*k4*trA3+2*k4*tlA3-6*tr*k4*tlA2+3*c45*trA2+3*c45*tlA2-6*tr*c45*tl+.. 

6*k4*tl*trA2)/hA2; 

k(5,5)=l/6*(2*k4*trA3-2*k4*tlA3+2*k5*trA3-2*k5*tlA3-3*c45*trA2-3*c45*tr2-3*c56*. 

trA2-3*c56*tlA2+6*tr*k4*tlA2+6*k5*tr*tlA2-6*j5*tr+6*j5*tl+6*tr*c45*tl+... 

6*tr*c56*tl-6*k4*tl*trA2-6*k5*trA2*tl)/hA2; 

k(5,6)=l/6*(-2*k5*trA3+2*k5*tlA3-6*k5*tt*tlA2+3*c56*trA2+3*c56*tlA2-6*tr*c56*tl+.. 

6*k5*trA2*tl)/hA2; 

k(5,10)=-l/6*(k4*trA3-k4*tlA3-3*k4*tl*trA2+3*c45*trA2+3*c45*tlA2+3*tr*k4*tlA2-6*... 

tr*c45*tl)/hA2; 

k(5,l  l)=-l/6*(6*tr*c45*tl-3*k5*tr*tlA2+3*k5*trA2*tl+6*tr*c56*tl-k5*trA3+k5*tlA3-... 
3*c56*trA2-3*c56*tlA2-3*tr*k4*tlA2+3*k4*tl*trA2-k4*trA3+k4*tlA3-3*c45*... 
trA2-3*c45*tlA2-6*j5*tr+6*j5*tl)/hA2; 

k(5,12)=-l/6*(k5*trA3-k5*tlA3+3*c56*trA2+3*c56*tlA2+3*k5*tr*tlA2-3*k5*trA2*tl-6*... 

tr*c56*tl)/hA2; 

k(6,5)=l/6*(-2*k5*trA3+2*k5*tlA3-6*k5*tr*tlA2+3*c56*trA2+3*c56*tlA2-6*tr*c56*tl+.. 

6*k5*trA2*tl)/hA2; 

k(6,6)=l/6*(2*k5*trA3-2*k5*tlA3+6*k5*tr*tlA2-3*c56*trA2-3*c56*tlA2-6*j6*tr+6*j6*... 

tl+6*tr*c56*tl-6*k5*trA2*tl)/hA2; 

k(6,ll)=-l/6*(k5*trA3-k5*tlA3+3*c56*trA2+3*c56*tlA2+3*k5*tr*tlA2-3*k5*trA2*tl-... 

6*tr*c56*tl)/hA2; 

k(6,12)=-l/6*(-k5*trA3+k5*tlA3-3*c56*trA2-3*c56*tlA2-3*k5*tr*tlA2+3*k5*trA2*tl-... 

6*j6*tr+6*j6*tl+6*tr*c56*tl)/hA2; 

k(7,l)=-l/6*(-kl*trA3+kl*tlA3+3*cl2*trA2+3*cl2*tlA2-3*kl*tr*tlA2+3*kl*trA2*tl-..-- 
6*j  1  *tr+6*j  1  *tl-6*c  1 2*tr*tl)/hA2; 

k(7,2)=l/6*(-kl*trA3+kl*tlA3+3*cl2*trA2+3*cl2*tlA2-3*kl*tr*tlA2+3*kl*trA2*tl-... 

6*cl2*tr*tl)/hA2; 

k(7,7)=l/6*(2*kl*trA3-2*kl*tlA3-6*kl*trA2*tl+3*cl2*trA2+3*cl2*tlA2-6*jl*tr+6*... 
j  1  *tl+6*kl  *tr*tlA2-6*c  12*tr*tl)/hA2; 

k(7,8)=-l/6*(2*kl*trA3-2*kl*tlA3-6*kl*trA2*tl+3*cl2*trA2+3*cl2*tlA2-6*cl2*tr*... 
tl+6*kl  *tr*tlA2)/hA2; 

k(8,l)=l/6*(-kl*trA3+kl*tlA3+3*cl2*trA2+3*cl2*tlA2-3*kl*tr*tlA2+3*kl*trA2*tl-... 

6*cl2*tr*tl)/hA2; 

k(8,2)=-l/6*(-3*k2*tr*tlA2+3*c2*trA2+3*c2*tlA2-6*j2rec(i)*tr+6*j2rec(i)*tl-6*j2*... 

tr+6*j2*tl-6*tl*tr*c23-6*tr*c2*tl+3*k2*tl*trA2-6*cl2*tr*tl+3*kl*trA2*... 

tI-k2*trA3+k2*tlA3+3*c23*trA2+3*c23*tlA2-kl*trA3+kl*tIA3+3*cl2*trA2+... 

3*c  1 2*tlA2-3*kl  *tr*tlA2)/hA2; 

k(8,3)=  1  /6*(-k2*trA3+k2*tlA3-3*k2*tr*tlA2+3  *c23  *trA2+3  *c23  *tlA2+3*k2*tl*trA2-... 
6*tl*tr*c23)/hA2; 

k(8,7)=-l/6*(2*kl*trA3-2*kl*tlA3-6*kl*trA2*tl+3*cl2*trA2+3*cl2*tlA2-6*cl2*tr*... 
tl+6*kl  *tr*tlA2)/hA2; 

k(8,8)=l/6*(6*k2*tr*tlA2+3*c2*trA2+3*c2*tlA2-6*j2rec(i)nr+6*j2rec(i)*tl-6*j2*... 
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tr+6*j2*tl-6*tl*tr*c23-6*tr*c2*tl-6*k2*tl*trA2-6*cl2*tr*tl-6*kl*trA2*... 
tl+2*k2*trA3-2*k2*tlA3+3*c23*trA2+3*c23*tlA2+2*kl*trA3-2*kl*tIA3+3*... 
c  12*trA2+3  *c  1 2*tlA2+6*k  1  *tr*tIA2)/hA2; 

k(8,9)=-l/6*(2*k2*trA3-2*k2*tlA3-6*k2*tl*trA2+3*c23*trA2+3*c23*tlA2-6*tl*tr*... 

c23+6*k2*tr*tlA2)/hA2; 

k(9,2)=l/6*(-k2*trA3+k2*tlA3-3*k2*tr*tlA2+3*c23*trA2+3*c23*tlA2+3*k2*tl*trA2-... 

6*tl*tr*c23)/hA2; 

k(9,3)=-l/6*(-3*k2*tr*tlA2-6*tl*tr*c23+3*k2*tl*trA2-3*k3*tr*tlA2+3*k3*trA2*tl-... 

6*tr*tl*c34-k2*trA3+k2*tlA3+3*c23*trA2+3*c23*tlA2-k3*trA3+k3*tlA3+... 

3*c34*trA2+3*c34*tIA2+3*c3*trA2+3*c3*tlA2-6*tl*c3*tr-6*j3*tr+6*j3*tI-... 

6*j3rec(i)*tr+6*j3rec(i)*tl)/hA2; 

k(9,4)=  1  /6*(-k3  *trA3+k3  *tlA3-3  *k3  *tr*tlA2+3*c34*trA2+3  *c34*tlA2+3  *k3*trA2*tl-. . . 
6*tr*tl*c34)/hA2; 

k(9,8)=-l/6*(2*k2*trA3-2*k2*tlA3-6*k2*tl*trA2+3*c23*trA2+3*c23*tlA2-6*tl*tr*... 
c23+6  *k2  *tr*tl  A2)/h  A2 ; 

k(9,9)=l/6*(6*k2*tr*tlA2-6*tl*tr*c23-6*k2*tl*trA2+6*k3*tr*tlA2-6*k3*trA2*tl-... 

6*tr*tI*c34+2*k2*trA3-2*k2*tlA3+3*c23*trA2+3*c23*tlA2+2*k3*trA3-... 

2*k3  *tlA3+3  *c34  *trA2+3*c34*tlA2+3  *c3  *trA2+3*c3  *tl  A2-6*tl*c3  *tr-. .. 
6*j3*tr+6*j3*tl-6*j3rec(i)*tr+6*j3rec(i)*tl)/hA2; 
k(9,10)=-l/6*(2*k3*trA3-2*k3*tlA3-6*k3*trA2*tl+3*c34*trA2+3*c34*tlA2-6*tr*... 
tI*c34+6*k3*tr*tlA2)/hA2; 

k(10,3)=l/6*(-k3*trA3+k3*tlA3-3*k3*tr*tlA2+3  *c34*trA2+3*c34*tlA2+3*k3*trA2*tl-... 
6*tr*tl*c34)/hA2; 

k(10,4)=-l/6*(-6*tr*c4*tl-6*tr*c45*tl-3*k3*tr*tlA2+3*k3*trA2*tI-6*tr*tl*c34-... 

k3*trA3+k3*tlA3+3*c34*trA2+3*c34*tlA2-3*tr*k4*tlA2+3*k4*tl*trA2-... 

k4*trA3+k4*tlA3+3*c45*trA2+3*c45*tlA2+3*c4*trA2+3*c4*tlA2-6*j4*tr+... 

6*j4*tl-6*j4rec(i)*tr+6*j4rec(i)*tl)/hA2; 

k(10,5)=l/6*(-k4*trA3+k4*tlA3+3*k4*tl*trA2+3*c45*trA2+3*c45*tlA2-3*tr*k4*tlA2-... 

6*tr*c45*tl)/hA2; 

k(10,9)=-l/6*(2*k3*trA3-2*k3*tlA3-6*k3*trA2*tl+3*c34*trA2+3*c34*tlA2-6*tr*tl*... 
c34+6*k3  *tr*tlA2)/h  A2 ; 

k(10,10)=l/6*(-6*tr*c4*tl-6*tr*c45*tl+6*k3*tr*tlA2-6*k3*trA2*tl-6*tr*tl*c34+2*... 

k3*trA3-2*k3*tlA3+3*c34*trA2+3*c34*tlA2+6*tr*k4*tlA2-6*k4*tI*trA2+... 

2*k4*trA3-2*k4*tIA3+3*c45*trA2+3*c45*tlA2+3*c4*trA2+3*c4*tlA2-6*j4*... 

tr+6*j4*tl-6*j4rec(i)*tr+6*j4rec(i)*tl)/hA2; 

k(10,ll)=-l/6*(2*k4*trA3-2*k4*tlA3-6*k4*tl*trA2+3*c45*trA2+3*c45*tlA2-6*tr*c45*... 

t!+6*tr*k4*tlA2)/hA2; 

k(l  l,4)=l/6*(-k4*trA3+k4*tlA3+3*k4*tl*trA2+3*c45*trA2+3*c45*tlA2-3*tr*k4*tlA2-... 
6*tr*c45*tl)/hA2; 

k(  1 1 ,5)=-  l/6*(-6*tr*c45*tl-3*k5  *tr*tl  A2+3*k5  *trA2*tl-6*tr*c56*tl-k5*trA3+k5  *tlA3+. . . 
3*c56*trA2+3*c56*tlA2-3*tr*k4*tlA2+3*k4*tl*trA2-k4*trA3+k4*tlA3+3*c45*... 
trA2+3  *c45  *tl  A2-6  *j5  *tr+6*j  5  *tl)/h  A2; 

k(ll,6)=l/6*(-k5*trA3+k5*tlA3+3*c56*trA2+3*c56*tlA2-3*k5*tr*tlA2+3*k5*trA2*tl-... 

6*tr*c56*tl)/hA2; 

k(ll,10)=-l/6*(2*k4*trA3-2*k4*tlA3-6*k4*tl*trA2+3*c45*trA2+3*c45*tlA2-6*tr*c45*tl+... 

6*tr*k4*tlA2)/hA2; 

k(l  1,1  l)=l/6*(2*k4*trA3-2*k4*tlA3+2*k5*trA3-2*k5*tlA3+3*c45*trA2+3*c45*tlA2+3*c56*... 
trA2+3*c56*tlA2-6*k4*tl*trA2-6*k5*trA2*tl-6*j5*tr+6*j5*tl+6*tr*k4*tlA2+... 
6*k5*tr*tlA2-6*tr*c45*tl-6*tr*c56*tl)/hA2; 

k(ll,12)=-l/6*(2*k5*trA3-2*k5*tlA3-6*k5*trA2*tl+3*c56*trA2+3*c56*tlA2-6*tr*c56*tl+... 

6*k5*tr*tlA2)/hA2; 

k(12,5)=l/6*(-k5*trA3+k5*tIA3+3*c56*trA2+3*c56*tlA2-3*k5*tr*tlA2+3*k5*trA2*tl-... 

6*tr*c56*tl)/hA2; 

k(12,6)=-l/6*(-k5*trA3+k5*tlA3+3*c56*trA2+3*c56*tlA2-3*k5*tr*tlA2+3*k5*trA2*tl-... 

6*j6*tr+6*j6*tl-6*tr*c56*tl)/hA2; 
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k(12,l  l)=-l/6*(2*k5*trA3-2*k5*tlA3-6*k5*trA2*tl+3*c56*trA2+3*c56*tIA2-6*tr*c56*tl+... 
6*k5*tr*tlA2)/hA2; 

k(12,12)=l/6*(2*k5*trA3-2*k5*tlA3-6*k5*trA2*tl+3*c56*trA2+3*c56*tlA2-6*j6*tr+6*j6*... 

tl+6*k5*tr*tlA2-6*tr*c56*tl)/hA2; 

Programs  used  to  calculate  cylinder  gas  torques,  given  angular  position  data 


at  0i  and  65. 


%  SPRESSF2 

% - - - 

%  Program  to  determine  cylinder  pressures  in  3  cylinder 

%  two  stroke  diesel  engine,  given  instantaneous  angular 

%  velocity  at  two  of  the  six  degrees  of  freedom 
%  (crankshaft  nose  and  flywheel) 


ic  =  [0  0  0  -1.2  -1.2  -10.4]*le-04;  %  set  initial  conditions 

% - 

%  load  data 


load  w51kl00.md 


t=[0;w51kl00(l:7938,l)]; 
tt  =  diff  (t); 

tt  =  reshape(tt,  126,63); 
mdt  =  mean(tt); 
for  i  =  1:63; 

tt(:,i)  =  tt(:,i)./mean(tt(:,i)); 


%  COMPUTES  THE  DELJT’S 
%  This  phase  locks  one  cycle 
%  Mean  dt  for  each  cycle 
%  Reference  each  cycle  to  its  own  mean 


end 

tt  =  tt**mean(mdt); 

dtrat  =  tt(l)/mean(tt(2:63,l)); 

teeth  =  linspace(0,2*pi,127);  tooth  =  2*pi/127; 

th5 1  =  ic(5)  +  tooth*dtrat; 

tttt  =  mean(tt);  %  This  ensemble  averages  the  phases 

cyctm  =  sum(tttt);  %  Time  for  one  complete  cycle 


thetarS  =  [ic(5)  teeth(l:125)+th51]; 

thetar5(127)  =  2*pi+ic(5); 

tjme5  =  [0,cumsum(tttt)-(l  -dtrat)  *tttt(l)]; 

time5(127)  =  cyctm; 

load  wllklOO.md 

t=  [0;  w  1 1  k 1 00(  1 :7920, 1 )] ; 

tt  =  diff(t); 

tt  =  reshape(tt,720,l  1)’; 
tt  =  tt(2:ll,:); 
tttt  =  mean(tt); 

thetarl  =  linspace(0,(2*pi),721); 
timel  =  [0  cumsum(tttt)]; 
timel  =  timel  *(cyctm/max(timel)); 
omega  =  (l/cyctm)*2*pi; 

N  =  512; 


%  Phase  lock  one  cycle 

%  Ensemble  average  the  phases 
%  Known  positions  of  the  O.E.  windows 

%  Adjust  times  to  agree 

%  set  number  of  nodes  for  solution 


%  calculate,  interpolate,  and  filter  position  data 

% - 

time  =  linspace(0,  cyctm,  N)’; 

dt  =  cyctm/N; 

theta  =  linspace(0,2*pi,N)’; 
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fcv  =  100*dt;  %  filtering  cut-off  value  as  ratio  of  sampling  freq 

thetarl  =  interpl  (time  1,  thetarl, time, ’spline5); 

thetar5  =  interpl  (time5,thetar5,  time,  ’spline7); 

thetal  =  theta  +  vfilt(thetarl -theta, fcv); 

theta5  =  theta  +  vfilt(thetar5-theta,fcv); 

% - 

%  derive  velocity  and  acceleration  data 

% . 

omegal  =  deriv(thetal,dt); 
omega5  =  deriv(theta5,dt); 
accell  =  deriv(omegal,dt); 
acceI5  =  deriv(omega5,dt); 

% - 

%  plot  raw  and  filtered  omegas  for  the  two  dofs 

% - 

figure  (1) 

plot  (time,  deriv(thetarl,dt),  ’mX),  grid,  hold  on 

plot  (time,  omegal,  V) 

plot  (time,  deriv(thetar5,dt),  ’cX5) 

plot  (time,  omega5,  ’b’) 

plot  (time,  ones(size(omega5))*omega,’g’) 

title  (Raw  and  Filtered  Angular  Velocity  at  DOFS  1  and  5*) 

xlabel  (’time  (seconds)5) 

ylabel  (’angular  velocity  (rad/sec)5) 

orient  landscape 

% . 

%  plot  raw  and  filtered  thetas  for  the  two  dofs 

% - 

figure  (2) 

plot  (time,  thetarl -theta,  ’mX),  grid,  hold  on 
plot  (time,  thetal -theta,  ’r7) 
plot  (time,  thetar5-theta,  ’cX) 
plot  (time,  theta5-theta,  ’b’) 

title  (Raw  and  Filtered  phase  deviation  at  DOFS  1  and  5*) 

xlabel  (’time  (seconds)5) 

ylabel  (’phase  deviation  from  mean  (rad)5) 

orient  landscape 

% - 


%  set  calibration  data 

C7„ 

j  1=0.02443; 

%lb*in*secA2/rad 

j2=0.2482; 

%lb*in*secA2/rad 

j3=0.1462; 

%lb*in*secA2/rad 

j4=0.2482; 

%lb*in*secA2/rad 

j5=7.222; 

%lb*in*secA2/rad 

j6=0.2870; 

%lb*in*secA2/rad 

k 1 =3 . 1 1 e6 ; 

%lb*in/rad 

k2=7.00e6‘, 

%lb*in/rad 

k3=7.00e6; 

%lb*in/rad 

k4=10.82e6; 

%lb*in/rad 

k5=1.304e6; 

%lb*in/rad 

cl2=0.01; 

%lb*in*sec/rad 

c23=0.01; 

%lb*in*sec/rad 

c34=0.01; 

%lb*in*sec/rad 

c45=0.01; 

%lb*in*sec/rad 
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c56=0.01; 
c2=0.013; 
c3=0.013; 
c4=0.013; 
Tload  =  1200; 
Fpar  =  100; 
Taux  = 168; 
Tpmp  =  9.5; 


%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%  lbf*in  Load  torque 

%  lbf  Parasitic  force 

%  lbf*in  Valvetrain  and  auxilliary  torque 

%  lbf*in  Oil  pump  torque 


% - 

%  Solve  equations  6,  5,  and  1 


global  AA  BB  CC  DD  TT 
icomega=  105.5*[1  11111]; 

thetar2=zeros(N,l);  thetar4=thetar2;  thetar6=thetar2;  thetar3=thetar2; 
omega2=zeros(N,l);  omegar3=omega2;  omega4=omega2;  omegar6=omega2; 
thetar3(l)  =  ic(3); 
omegar3(l)  =  icomega(3); 

% - solve  equation  6  for  theta6  and  omega6 - 

bcval  =  [ic(6),ic(6)+2*pi]; 

AA=j6;  BB-c56;  CC  =  k5; 

DD  =  -Tload  +  c56*omega5  +k5*theta5; 
ff  =  zeros(N,l); 
kk  =  zeros(N,N); 
for  i  =  1:(N-1); 
h  =  time(i+l)-time(i); 

k  =  -(AA/h)*[l  -1;-1  1]  +  (BB/2)*[-l  1;-1  1]  +  (CC*h/6)*[2  1;1  2]; 
f=(h/2)*[DD(i);DD(i+l)]; 
kk(i:i+l,i:i+l)  =  kk(i:i+l,i:i+l)  +  k; 
ff(i:i+l)  =  ff(i:i+l)  +  f; 
end 

% - apply  boundary  conditions - 

kk(l,:)  =  zeros(l,N);  kk(N,:)  =  zeros(l,N); 

kk(  1,1)  =  1;  kk(N,N)=  1; 

ff(l)  =  bcval(l);  ff(N)  =  bcval(2); 

% - solve  matrix  eqn - 

theta6  =  kk\ff; 

omega6  =  deriv(theta6,dt); 

% - solve  equation  5  for  theta4  — 

thetar4  =  (Taux  +  k5*(theta5-theta6)  +  c56*(omega5-omega6)  +... 

k4*theta5  +  J5*accel5)./k4; 
theta4  =  theta  +  vfilt(thetar4-theta,fcv); 

% . -  solve  equation  1  for  theta2 - 

thetar2  =  (Tpmp  +  kl*thetarl  +  jl*accell)./kl; 
theta2  =  theta  +  vfilt(thetar2-theta,fcv); 

% - compute  omega  and  accel  for  dofs  2  and  4 - 

omega2  =  deriv(theta2,dt);  %  raw  velocities 

omega4  -  deriv(theta4,dt); 

accel2  =  deriv(omega2,dt);  %  calculated  accelerations 

acce!4  =  deriv(omega4,dt); 


%  Solve  equations  2,  3,  and  4  in  three  steps 

% - - 

R  =  2.25;  %  Crankshaft  eccentricity  (in) 

W  =  7.556;  %  reciprocating  weight(lbf) 

g  =  386;  %  The  acceleration  of  gravity  (in/secA2) 
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B  =  3.875;  %  Cylinder  Bore  (in) 

L  =  8.80;  %  Connecting  Rod  length  (in) 

j2rec=(W*RA2/(2*g))*(l-cos(2*theta));  %lb*in*secA2 

j3rec^W*RA2/(2*g))*(l-cos(2*(theta-4*pi/3)));  %lb*in*secA2 

j4rec=(W*RA2/(2*g))*(l-cos(2*(theta>2*pi/3)));  %Ib*in*secA2 

Tlcyl  =  zeros(N,l);  T2cyl=Tlcyl;  T3cyl  =  Tlcyl; 
pref  =  14.706  +  0.00205*((omega*60/(2*pi))-634);  %  ref  press  (psia) 

N1  =  min(find(thetal>=(2*pi/3)));  %  index  for  TDC  Cyl  #3 

N2  =  min(find(thetal>=(4*pi/3)));  %  index  for  TDC  Cyl  #2 

accel3  =  zeros(size(accel2»; 

si  =  R*cos(theta)  +  sqrt(LA2  -  (RA2)*sin(theta).A2); 

s2  =  R*cos(theta-4*pi/3)  +  sqrt(LA2  -  (RA2)*sin(theta~4*pi/3).A2); 

s3  =  R*cos(theta-2*pi/3)  +  sqrt(LA2  -  (RA2)*sin(theta-2*pi/3)  A2); 

Spl  =  deriv(sl,dt);  Spdl  =  deriv(Spl,dt);  %  Piston  speed  (in/sec)  and 

Sp2  =  deriv(s2,dt);  Spd2  =  deriv(Sp2,dt);  %  Piston  acceleration  (in/secA2) 

Sp3  =  deriv(s3,dt);  Spd3  =  deriv(Sp3,dt); 

Tlrec  =  -(W/g)*Spdl.*(Spl/omega);  %  in*lbf  Reciprocating  torque 

T2rec  =  -(W/g)*Spd2.*(Sp2/omega); 

T3rec  =  -(W/g)*Spd3.*(Sp3/omega); 

Tparl  =  Fpar*abs(Spl/omega);  %  in*lbfParasitic  torque 

Tpar2  =  Fpar*abs(Sp2/omega); 

Tpar3  =  Fpar*abs(Sp3/omega); 

% - step  one:  Determine  known  values  of  Tcyl  from  pref - 

% - (known  values  of  Tcyl  are  0) - 

% - Step  two:  Solve  for  theta3  throughout  cycle - 

% - solve  equation  2  for  theta3 - 

thetar3(Nl:(N2-l))  =  ((j2+j2rec(Nl:(N2-l)))  *accel2(Nl:(N2-l))  +  ... 
c  1 2*(omega2(N  1  :(N2- 1  ))-omegal  (N 1  :(N2- 1 )))  +  ... 
k  1  *(theta2(N  1 :  (N2- 1  ))-theta  1  (N 1  :(N2- 1 )))  +  k2*theta2(Nl:(N2-l))  +  ... 
c2*omega2(N  1  :(N2- 1 ))  -  Tlcyl(Nl:(N2-l))  -  Tlrec(Nl:(N2-l))  +  ... 

Tparl(Nl:(N2-l)))./k2; 

% - solve  equation  3  for  theta3  and  omega3 - 

BB  =  c23+c34+c3;  CC  =  k2+k3; 

DDT  =  T2cyl  +  T2rec  -  Tpar2  +  c23*omega2  +  k2*theta2  +  c34*omega4  +  k3*theta4; 
for  i  =  l:(Nl-2); 

AA  =  j3  +  j3rec(i:i+l);  TT  =  time(i:i+l);  DD  =  DDT(i:i+l); 

[T,X]  =  ode45(’seqns2,,[time(i),time(i+l)],[thetar3(i);omegar3(i)]); 
thetar3(i+l)  =  X(length(T),l); 
omegar3(i+l)  =  X(length(T),2); 
end 

% - solve  equation  4  for  theta3 - 

thetar3(N2:N)  =  ((j4+j4rec(N2:N))*accel4(N2:N)  +  k3*theta4(N2:N)  +  ... 

(c45+c4)*omega4(N2:N)  -  c45*omega5(N2:N)  +  k4*(theta4(N2:N)  - ... 
theta5(N2:N))  -  T3cyl(N2:N)  -  T3rec(N2:N)  +  Tpar3(N2:N))./k3; 

% - filter  theta3  and  derive  omega3  and  accel3 - 

theta3  =  theta  +  vfilt(thetar3-theta,fcv); 
omega3  =  deriv(theta3,dt); 
accel3  =  deriv(omega3,dt); 

% . Step  three:  solve  for  remaining  Tcyl  values - 

% - solve  equation  2  for  Tlcyl - 

Tlcyl(l:(Nl-l))  =  -Tlrec(l:(Nl-l))  +  Tparl(l:(NM))  +  (j2+j2rec(l:(NM))).*acceI2(l:(Nl-l))  +  . 
cl2*(omega2(l:(Nl-l))-omegal(l:(Nl-l)))  +  kl*(theta2(l:(Nl-l))-thetal(l:(Nl-l»)  +... 
c23*(omega2(  1  :(N  1  - 1  ))-omega3(  1  :(N1  - 1 )))  +  k2*(theta2(l:(NM))-theta3(l:(Nl-l)))  +... 
c2*omega2(  1  :(N1  - 1 )); 

Tlcyl(N2:N)  =  -Tlrec(N2:N)  +  Tparl(N2:N)  +  (j2+j2rec(N2:N))  *acce!2(N2:N)  +  ... 
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cl2*(omega2(N2:N)-omegal(N2:N))  +  kl*(theta2(N2:N)-thetal(N2:N))  +  ... 
c23*(omega2(N2:N)-omega3(N2:N))  +  k2*(theta2(N2:N)-theta3(N2:N))  +  ... 
c2*omega2(N2:N); 

% . -  solve  equation  3  for  T2cyl - 

T2cyl(Nl:N)  =  -T2rec(Nl:N)  +  Tpar2(Nl:N)  +  (j3+j3rec(Nl:N)).*accel3(Nl:N)  +  ... 
c23*(omega3(Nl:N)-omega2(Nl:N))  +  k2*(theta3(Nl:N)-theta2(Nl:N))  +  ... 
c34*(omega3(Nl:N)-omega4(Nl:N))  +  k3*(theta3(Nl:N)-theta4(Nl:N))  +  ... 
c3*omega3(N  1  :N); 

% - solve  equation  4  for  T3cyl - 

T3cyl(l:(N2-l))  =  -T3rec(l:(N2-l))  +  Tpar3(l:(N2-l))  +  (j4+j4rec(l:(N2-l)))  *acce!4(l:(N2-l))  +  . 
c34*(omega4(  1  :(N2- 1  ))-omega3(  1  :(N2- 1 )))  +  k3*(theta4(l:(N2-l))-theta3(l:(N2-l)))  +... 
c45  *(omega4(  1  :(N2- 1  ))-omega5  ( 1  :(N2- 1 )))  +  k4*(theta4(l:(N2-l))-theta5(l:(N2-l)))  +... 
c4*omega4(  1  :(N2- 1 )); 

% - Convert  to  FT*LBF - 

Tlcyl  =  Tlcyl./12; 

T2cyl  =  T2cyl./12; 

T3cyl  =  T3cyl./12; 

% - 

%  Determine  and  plot  measured/predicted  cylinder  torques 

% - 

load  walklOO.md  %  ECA  cylinder  #1  pressure  data 

load  wb  1  kl  OO.md  %  ECA  cylinder  #2  pressure  data 

load  wclkl OO.md  %  ECA  cylinder  #3  pressure  data 

pa  =  reshape  (walkl00,5,720);  Plcyl  =  [pa(l,:)  pa(l)]; 
pb  =  reshape  (wblkl 00,5, 720);  P2cyl  =  [pb(l,:)  pb(l)]; 
pc  =  reshape  (wclkl00,5,720);  P3cyl  =  [pc(l ,:)  pc(l)]; 
pos  =  linspace(0,2*pi,721); 
dtm  =  (2*pi/omega)/720; 

slm  =  R*cos(pos)  +  sqrt(LA2  -  (RA2)*sin(pos).A2); 

s2m  =  R*cos(pos-4*pi/3)  +  sqrt(LA2  -  (RA2)*sin(pos-4*pi/3)  A2); 

s3m  =  R*cos(pos-2*pi/3)  +  sqrt(LA2  -  (RA2)*sin(pos-2*pi/3)  A2); 

Splm  =  deriv(slm,dtm);  %  Piston  speed  (in/sec) 

Sp2m  =  deriv(s2m,dtm); 

Sp3m  =  deriv(s3m,dtm); 

Tlcylm  = -((Plcyl. *pref-pref)*(pi*BA2/4).*Splm/omega)/12;  %  (ft*lbf) 

T2cylm  =  -((P2cyl.*pref-pref)*(pi*BA2/4).*Sp2m/omega)/12;  %  (ft*lbf) 

T3cylm  =  -((P3cyl.*pref-pref)*(pi*BA2/4).*Sp3m/omega)/12;  %  (ft*Ibf) 
figure(3) 

plot(pos,Tlcylm,,bX,,pos,T2cylm,,bO,,pos,T3cylm,’b.,) 
legend(’Cyl  #1  VCyl  #2’,Cyl  #30 
grid,  hold  on 

plot(theta  1  ,T  1  cy  1 ,  ’rX  ’,  theta  1  ,T2cy  1 ,  VO  ’,  theta  1  ,T3  cy  1 ,  ’r .  0 

title( Measured  (blue)  and  Predicted  (red)  Cylinder  Gas  Torques’) 

ylabel(Torque  (ft*lbf)0 

xlabel(’Crank  Angle*) 

orient  landscape 

figure  (4) 

plot(pos,T  1  cylm+T2cylm+T3cy  lm,  b*) 
grid,  hold  on 

plot(thetal  ,T1  cyl+T2cyl+T3cyl,  V*) 

title( Measured  (blue)  and  Predicted  (red)  Total  Cylinder  Gas  TorqueO 

ylabel  (Torque  (ft*lbf)0 

xlabel  (’Crank  AngleO 

orient  landscape 

figure(5) 
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degm=  180*pos/pi; 

degp  =  180*theta/pi; 

subplot(3,l,l) 

plot(degm,T  1  cy  lm,  Tc.  ) 

grid,  hold  on 

plot(degp,Tlcyl,’k) 

legend  (Meas  ’,Pred  ) 

ylabel  (’Cyl  #1  Torque  (ft*lbf)) 

title(’1000  RPM,  100  ft*lbf  Cylinder  Gas  Torques) 

axis([0, 360, -500, 1000]) 

subplot(3,l,2) 

plot(degm,T2cylm,  k. ) 

grid,  hold  on 

plot(degp,T2cyl,k) 

ylabel  (Gy  1  #2  Torque  (ft*lbf)) 

axis([0, 360, -500, 1000]) 

subplot(3,l,3) 

plot(degm,T3cylm,  k. ) 

grid,  hold  on 

plot(degp,T3cyl,k) 

ylabel  (’Cyl#3  Torque  (ft*lbf)) 

xlabel(’Crank  Angle  (deg)) 

axis([0,360,-500, 1 000]) 

orient  tall 

figure  (6) 

plot(degm,T  1  cy  lm+T2cyim+T3cylm,  k. ) 
grid,  hold  on 

plot(degp,Tlcyl+T2cyl+T3cyl,k) 

title  ('1 000  RPM,  100  ft*lbf  Total  Cylinder  Gas  Torque) 

legend  (Meas  ’,Pred  ) 

ylabel  (Total  Cylinder  Gas  Torque  (ft*lbf)) 

xlabel  (’Crank  Angle  (deg)) 

axis([0,360,-400,800]) 

orient  tall 

%  DERIV 

%  Function  to  determine  1-D  derivative  of  a  vector  using 

%  a  central  difference  technique. 

%  Xd  =  Deriv(X,dt) 

%  Returns  the  derivative  of  the  vector  X  as  a  function  of 

%  t,  given  the  time  step  dt.  Default  value  for  dt  is  1 . 

function  xd  =  deriv(x,dt) 
if  nargin  ==  1 ; 
dt  =  1; 
end 

N  =  length(x); 
xd  =  zeros(size(x)); 
xdf  =  diff(x); 

xd(2:(N-l))  =  (xdf(2:(N-l))  +  xdf(l:(N-2)))./(2*dt); 
xd(l)  =  xdf(l)/dt; 
xd(N)  =  xdf(N-l)/dt; 

%  Forward  difference  format 
%xd(l:N-l)  =  xdf/dt; 

%xd(N)=xd(l); 
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%  SEQNS2 

%  Function  to  solve  arbitrary  2nd  order  ode  in  form 

%  AAxdd  +  BBxd  +  CCx  =  DD 

%  where  A  A,  BB,  CC,  and  DD  are  global  variables 

function  xdot=seqns2(t,x) 

global  AA  BB  CC  DD  TT 

DDS  =  DD(1)  +  d  iff(DD)  *(  t-TT  ( 1  ))/d  i  ff(TT) ; 

AAS  =  AA(1)  +  diff(AA)*(t-TT(l))/diff(TT); 
xdot  =  zeros(2,l); 
xdot(l)  =  x(2); 

xdot(2)  =  (DDS  -  BB.*x(2)  -  CC.*x(l))./AAS; 

%  VFILT 

%  Function  performs  fast  fourier  transform  filtration 

%  of  high  frequency  components  of  given  data 

%  Y  =  VFILT(X,FCV) 

%  Y  is  filtered  data 

%  X  is  input  data 

%  FCV  is  frequency  cutoff  value  (as  a  fraction  of 
%  the  sampling  frequency) 


function  Y  =  vfilt(x,fcv) 

N  =  length(x); 

%  number  of  data  points 

D  =  fft(x); 

%  Fourier  transform  of  data 

ico  =  fix(fcv*N); 

%  index  in  D 

DF  =  zeros(l,N); 

DF(l:(ico+l))  =  D(l:(ico+l)); 

%  filter  high  frequencies 

DF((N-ico+l):N)  =  D((N-ico+l):N); 

%  mirror  image  data 

Y  =  ifft(DF); 

%  inverse  fft 

Y  =  real(Y); 

Y  =  reshape(Y,size(x)); 

Program  used  for  analysis 

of  torsional  model  natural  frequencies  and  modal 

shapes. 

% - 

%  NATFREQ 


Determines  natural  frequencies  of  torsional  model 


j  1=0.02443; 

J2=0.2482; 

j3=0.1462; 

j4=0.2482; 

j5=7.222; 

j6=0.2870; 

jrec=0.02478; 

kl=3.11e6; 

k2=7.00e6; 

k3=7.00e6; 

k4=10.82e6; 

k5=1.304e6; 

cl  2=0.01; 

c23=0.01; 

c34=0.01; 

c45=0.01; 


%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in*secA2/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 

%lb*in*sec/rad 
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c56=0.01;  %lb*in*sec/rad 

c2=0.013;  %Ib*in*sec/rad 

c3=0.013;  %lb*in*sec/rad 

c4=0.013;  %lb*in*sec/rad 

% - 

J  =  [jl  0  0  0  0  0;0  j2+jrec  0  0  0  0;0  0j3+jrec  0  0  0;... 

0  0  0  j4+jrec  0  0;0  0  0  0  j5  0;0  0  0  0  0  j6]; 

% . 

K  =  [kl  -kl  0  0  0  0;-kl  kl+k2  -k2  0  0  0;0  -k2  k2+k3  -k3  0  0;... 
0  0  -k3  k3+k4  -k4  0;0  0  0  -k4  k4+k5  -k5;0  0  0  0  -k5  k5]; 

% - - - 

ws  =  eig(K/J); 

w  =  sqrt(ws)  %  (rad/sec)  natural  frequencies 

whz  =  w/(2*pi)  %  (Hz)  natural  frequencies 

for  i  =  1:6; 
kj  =  K-J*w(i)A2; 

A(l,i)  =  1; 

A(2,i)  =  -A(l,i)*kj(l,l)/kj(l,2); 

A(3,i)  =  (-A(  1  ,i)*kj(2, 1  )-A(2,i)*kj(2,2»/kj(2,3); 

A(4,i)  =  (-A(2,i)*kj(3,2)-A(3,i)*kj(3,3»/kj(3>4); 

A(5,i)  =  (-A(3,i)*kj(4,3)-A(4,i)*kj(4,4))/kj(4,5); 

A(6,i)  =  (-A(4,i)*kj(5,4)-A(5,i)*kj(5,5))/kj(5,6); 
end 

[whz,I]  =  sort(abs(whz)); 

A 

figure  (1) 
for  i  =  2:6; 

subplot(5,l,i’l) 
plot(A(:,I(i)),’k’),  grid,  hold  on 
plot(A(:J(i)),’ko’) 
ylabel  ([num2str(whz(i),4),’Hz’]) 

%title  ([Mode  ’,num2str(i,l)]) 
end 

orient  tall 
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APPENDIX  F.  UNCERTAINTY  ANALYSIS 


Various  sources  of  error  will  affect  the  accuracy  of  results  in  this  study.  The 
following  discussion  contains  a  qualitative  summary  of  the  potential  sources  of 
experimental  uncertainty. 

1.  Optical  Encoder 

The  use  of  the  flexible  coupling  results  in  a  torsional  oscillation  of  the  encoder 
disk  about  the  actual  angular  position  of  the  crankshaft  nose  (0i).  The  frequency  of  the 
oscillation  is  apparent  from  the  measured  data  and  discussed  in  Appendix  C.  Heidenhain 
[Ref  21]  lists  a  kinematic  error  of  transfer  of  ±10”  for  the  encoder  coupling, 
corresponding  to  a  vibrational  amplitude  of  about  1  x  10'4  radians.  The  oscillation  seen 
in  experimental  data  varies  from  about  that  value  up  to  1  x  10'3  radians.  But  this  high 
frequency  oscillation  is  easily  filtered  out  of  the  raw  data,  and  its  amplitude  is  still  an 
order  of  magnitude  smaller  than  the  amplitude  of  the  crankshaft  angular  velocity 
fluctuations. 

The  3600  count  optical  encoder  allows  a  theoretical  adjustment  to  within  0.1°,  but 
this  is  only  as  accurate  as  the  TDC  alignment  for  the  engine.  Due  to  an  inadequate 
method  of  determining  TDC,  this  error  may  grow  to  1  or  2°. 

2.  Flywheel 

The  TDC  position  for  the  flywheel  data  is  determined  by  a  comparison  of  the  first 
At  to  subsequent  values,  assuming  that  the  crankshaft  has  zero  twist  at  TDC.  This  is  a 
reasonable  assumption,  but  not  completely  accurate.  Therefore,  an  expected  error  is 
introduced  due  to  ambiguity  in  determination  of  the  flywheel  TDC  angular  position. 
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Tooth  to  tooth  variation  will  result  in  a  cyclic  error  for  the  flywheel  data.  A 
motoring  measurement  of  the  flywheel  at  a  constant  speed  could  be  used  to  generate  a 
correction  signal,  eliminating  this  error  [Ref  8].  This  was  not  done  for  this  study, 
assuming  the  error  would  be  small  relative  to  the  measured  velocity  fluctuations. 

A  radial  vibration  is  present  in  the  flywheel  during  engine  operation.  Since  the 
proximeter  is  mounted  to  view  the  ring  gear  teeth  from  the  side,  this  radial  vibration  will 
cause  the  position  of  the  proximeter  relative  to  the  tooth  to  oscillate  up  and  down.  This 
up  and  down  oscillation  will  induce  a  cyclic  variation  in  the  pulse  width  not  associated 
with  crankshaft  rotational  velocity  fluctuations.  However,  it  is  expected  that  this 
variation  is  reasonably  small  enough  to  be  ignored. 

3.  Measurement  Error 

Engine  speed  and  load  vary  slowly  during  data  collection.  Since  data  for  the 
various  runs  are  not  collected  simultaneously,  there  are  small  differences  in  speed  and 
load  for  data  comprising  a  set.  Typically,  engine  speed  variation  was  within  20  RPM  of 
the  stated  value,  and  load  was  within  1  ft*lbf.  A  correction  is  made  in  the  pressure 
prediction  program  to  account  for  small  differences  in  rotation  speed  between  the  0]  and 
05  data. 
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